diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index 6912350..de90298 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -39,9 +39,9 @@ is committed and listed here with its exports. | assets/icons/interaction.contact.general.svg, assets/icons/interaction.contact.general.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.contact.general.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/interaction.contact.pair.svg, assets/icons/interaction.contact.pair.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.contact.pair.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/interaction.fsi.svg, assets/icons/interaction.fsi.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.fsi.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | -| assets/icons/interaction.mass.svg, assets/icons/interaction.mass.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.mass.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | +| assets/icons/model.mass.create.svg, assets/icons/model.mass.create.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.mass.create.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/interaction.rigid_wall.svg, assets/icons/interaction.rigid_wall.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.rigid_wall.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | -| assets/icons/interaction.spring.svg, assets/icons/interaction.spring.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.spring.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | +| assets/icons/model.spring.create.svg, assets/icons/model.spring.create.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.spring.create.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/interaction.tie.svg, assets/icons/interaction.tie.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.tie.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/mesh.delaunay.points.svg, assets/icons/mesh.delaunay.points.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.delaunay.points.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/mesh.element.add.svg, assets/icons/mesh.element.add.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.element.add.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | diff --git a/ROADMAP.md b/ROADMAP.md index 4cb66c1..da4cbb7 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -22,10 +22,10 @@ hex meshes today and it gives the best elements where stresses are read. | M1 Geometry in and surface preparation | Open `.solidcraft` designs with persistent face names, assemblies placed (done); face-tagged 3MF (done); test corpus, `cargo xtask corpus` (done); exact predicates (done); quality metrics and certified validity gate (done); weld, closed/manifold/orientation and self-intersection checks, feature edges (`modelcraft-mesh::surface`, `mesh.surface.check`, done) | in progress | | M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity; the `modelcraft_mesh::size::SizeField` trait and edge-division helpers are in); Bowyer–Watson Delaunay (exact, symbolic perturbation, verified; `modelcraft-mesh::delaunay`, done); robust CDT boundary recovery (Diazzi et al. 2023: implicit LNC Steiner points on segments, gift-wrapping facet recovery, inside/outside labelling, B-rep face tags on boundary tet faces, few-ulp rounding repair checked by the validity gate; `modelcraft-mesh::cdt` and `mesh.tet.generate`, done: all 20 valid corpus parts); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress | | M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input (done: `deck.write`), Nastran SOL 101/103 (done: `deck.write`) and LS-DYNA implicit (done); the `.modelcraft` file (text definition, binary mesh, geometry references reopened on load, journal; atomic saves, autosave and recovery; undo/redo on the journal: done) | in progress | -| M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion; 2D elements: plane stress, plane strain and axisymmetric tri3, tri6, quad4 (incompatible modes QM6), quad8 with two degrees of freedom per node, centrifugal loads, `mesh.planar` 2D meshes of polygons and planar faces, verified by patch tests, Kirsch, Lamé (plane strain and axisymmetric), rotating discs and cylinders, and a plane stress cantilever; deck mapping to Abaqus CPS/CPE/CAX, Nastran membranes, LS-DYNA ELFORM 12/13); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started | +| M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion; 1D elements I: rod2/rod3 trusses, point masses, rotary inertia, springs and dashpots in both solvers and all three decks, verified by the method of joints, the three-bar truss, spring–mass, two-mass chain and torsional pendulum frequencies, and damped decay; 2D elements: plane stress, plane strain and axisymmetric tri3, tri6, quad4 (incompatible modes QM6), quad8 with two degrees of freedom per node, centrifugal loads, `mesh.planar` 2D meshes of polygons and planar faces, verified by patch tests, Kirsch, Lamé (plane strain and axisymmetric), rotating discs and cylinders, and a plane stress cantilever; deck mapping to Abaqus CPS/CPE/CAX, Nastran membranes, LS-DYNA ELFORM 12/13); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started | | M5 Structured hex | 2D quad meshing of faces (`modelcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | | M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned | -| M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, springs, point masses; amplitudes, multiple steps; nonlinear static steps in decks | planned | +| M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, connectors; amplitudes, multiple steps; nonlinear static steps in decks | planned | | M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; heat transfer | planned | | M9 Post-processing and deck readers | Results readers for public formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`); XY and path plots, cuts, display groups, free-body forces; deck readers that keep unknown keywords; mesh convergence studies | planned | | M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `modelcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (X2: done — updated Lagrangian elements, Jaumann/Hughes–Winget and Green–Naghdi/Flanagan–Taylor rates, J2 bilinear/table plasticity, Johnson–Cook strength with heating and damage, erosion; materials gain `plastic`, `damage`, `erosion` with `model.material.plasticity`/`damage`; Dienes shear, rigid rotation, J2 uniaxial and Taylor-bar tests); contact (global search, penalty with friction, rigid walls); runs on hex and tet meshes | started | diff --git a/assets/icons/interaction.mass.vectorcraft b/assets/icons/interaction.mass.vectorcraft deleted file mode 100644 index 47cfda1..0000000 Binary files a/assets/icons/interaction.mass.vectorcraft and /dev/null differ diff --git a/assets/icons/interaction.spring.vectorcraft b/assets/icons/interaction.spring.vectorcraft deleted file mode 100644 index caf2653..0000000 Binary files a/assets/icons/interaction.spring.vectorcraft and /dev/null differ diff --git a/assets/icons/interaction.mass.svg b/assets/icons/model.mass.create.svg similarity index 95% rename from assets/icons/interaction.mass.svg rename to assets/icons/model.mass.create.svg index b4111b2..70f253a 100644 --- a/assets/icons/interaction.mass.svg +++ b/assets/icons/model.mass.create.svg @@ -1,6 +1,6 @@ - interaction.mass + model.mass.create diff --git a/assets/icons/model.mass.create.vectorcraft b/assets/icons/model.mass.create.vectorcraft new file mode 100644 index 0000000..c92d816 Binary files /dev/null and b/assets/icons/model.mass.create.vectorcraft differ diff --git a/assets/icons/interaction.spring.svg b/assets/icons/model.spring.create.svg similarity index 95% rename from assets/icons/interaction.spring.svg rename to assets/icons/model.spring.create.svg index 8e42a97..cadb03a 100644 --- a/assets/icons/interaction.spring.svg +++ b/assets/icons/model.spring.create.svg @@ -1,6 +1,6 @@ - interaction.spring + model.spring.create diff --git a/assets/icons/model.spring.create.vectorcraft b/assets/icons/model.spring.create.vectorcraft new file mode 100644 index 0000000..82a2de2 Binary files /dev/null and b/assets/icons/model.spring.create.vectorcraft differ diff --git a/assets/icons/src/interaction.mass.json b/assets/icons/src/model.mass.create.json similarity index 100% rename from assets/icons/src/interaction.mass.json rename to assets/icons/src/model.mass.create.json diff --git a/assets/icons/src/interaction.spring.json b/assets/icons/src/model.spring.create.json similarity index 100% rename from assets/icons/src/interaction.spring.json rename to assets/icons/src/model.spring.create.json diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index b1c10fd..03710c9 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -63,6 +63,18 @@ //! | tri3, quad4 | | S3 / S3R, S4 / S4R | | | //! | quad8 | | S8R (reduced only) | | | //! | beam2, rod2 | | | B31 | T3D2 | +//! | rod3 | | | | T3D3 | +//! +//! Trusses in a 2D model are T2D2 and T2D3. +//! +//! Point masses, springs and dashpots are elements numbered after the mesh's, one per resolved +//! copy, in an element set per definition (`-mass`, `-inertia`, `-spring`, made +//! unique against the model's sets): `MASS` with `*MASS` (one mass for every direction), +//! `ROTARYI` with `*ROTARY INERTIA` (I11, I22, I33, I12, I13, I23), and `SPRING1`/`DASHPOT1` to +//! the ground and `SPRING2`/`DASHPOT2` between two nodes along a global degree of freedom, or +//! `SPRINGA`/`DASHPOTA` along the line of the nodes, with `*SPRING`/`*DASHPOT`. A mass that +//! differs by direction and a spring along a skew fixed vector are reported (they need a local +//! orientation, not written yet). //! //! | ModelCraft | Plane stress | Plane strain | Axisymmetric | //! |---|---|---|---| @@ -137,13 +149,16 @@ pub fn write_resolved(model: &Model, r: &Resolved) -> Result { crate::refuse_fluid(FORMAT, model)?; let by_type = elements_by_type(model, r)?; check_expressible(model)?; + let discrete = discrete_groups(model, r)?; let mut d = Deck::default(); d.heading(model); d.nodes(model)?; d.elements(&by_type); + d.discrete_elements(&discrete); d.sets(model, r)?; d.materials(model)?; d.sections(model)?; + d.discrete_properties(&discrete)?; d.couplings(model, r); d.steps(model, r)?; Ok(d.out) @@ -198,6 +213,7 @@ pub fn element_type(kind: ElementKind, section: &Section) -> std::result::Result (S::Shell { .. }, K::Tri6, ..) => return Err("tri6 shells (Abaqus has no general-purpose six-node triangular shell)".into()), (S::Beam { .. }, K::Beam2, ..) => "B31", (S::Rod { .. }, K::Rod2, ..) => "T3D2", + (S::Rod { .. }, K::Rod3, ..) => "T3D3", (s, ..) => return Err(format!("{} elements in a {} section", kind.name(), s.name())), }; Ok(if f.hybrid { format!("{plain}H") } else { plain.to_string() }) @@ -216,7 +232,9 @@ fn elements_by_type<'m>(model: &'m Model, r: &Resolved) -> Result { count += elements.len(); by_type.entry(t).or_default().extend(elements); @@ -284,6 +302,90 @@ fn check_expressible(model: &Model) -> Result<()> { unsupported(FORMAT, bad) } +/// The property of a group of discrete elements. +enum DiscreteProperty { + Mass(f64), + Inertia([f64; 6]), + /// `*SPRING` or `*DASHPOT`: the degrees of freedom line (empty for the axial form) and the + /// stiffness or damping coefficient. + Spring { + dashpot: bool, + dofs: String, + value: f64, + }, +} + +/// The discrete elements of one definition: element type, element set, members and property. +struct DiscreteGroup { + element_type: &'static str, + elset: String, + members: Vec<(u32, Vec)>, + property: DiscreteProperty, +} + +/// Point masses, springs and dashpots as element groups (one per definition and element type, +/// labels counting up after the mesh's), or everything Abaqus cannot take. +fn discrete_groups<'r>(model: &Model, r: &'r Resolved) -> Result> { + let mut next = model.elements().keys().next_back().map_or(1, |l| l.saturating_add(1)); + let mut label = || -> Result { + let l = next; + if l > modelcraft_model::MAX_LABEL { + return Err(Error::Limit(format!( + "no element label left below {} for point masses, springs and dashpots", + modelcraft_model::MAX_LABEL as u64 + 1 + ))); + } + next += 1; + Ok(l) + }; + let mut taken: BTreeSet = model.sets().keys().map(|n| n.to_ascii_uppercase()).collect(); + let mut groups: Vec = Vec::new(); + // The (definition, element type) the last group holds. + let mut current: Option<(&str, &'static str)> = None; + let mut push = |name: &'r str, suffix: &str, element_type: &'static str, property: DiscreteProperty, member: (u32, Vec)| { + let key = (name, element_type); + if current.as_ref() == Some(&key) + && let Some(g) = groups.last_mut() + { + g.members.push(member); + return; + } + let elset = crate::unique_name(&format!("{name}-{suffix}"), &mut taken); + groups.push(DiscreteGroup { element_type, elset, members: vec![member], property }); + current = Some(key); + }; + let mut bad = BTreeSet::new(); + for m in &r.masses { + if m.mass.iter().any(|c| *c != 0.0) { + if m.mass[0] != m.mass[1] || m.mass[0] != m.mass[2] { + bad.insert(format!("point mass `{}`: a mass that differs by direction (needs a local orientation, not written yet)", m.name)); + } + push(&m.name, "mass", "MASS", DiscreteProperty::Mass(m.mass[0]), (label()?, vec![m.node])); + } + } + for m in &r.masses { + if m.inertia.iter().any(|c| *c != 0.0) { + push(&m.name, "inertia", "ROTARYI", DiscreteProperty::Inertia(m.inertia), (label()?, vec![m.node])); + } + } + for s in &r.springs { + let what = s.kind(); + let line = matches!(s.direction, modelcraft_model::ResolvedDirection::Line(_)); + let (element_type, dofs, nodes) = match (s.second, crate::axis_dof(s.direction), line) { + (Some(b), _, true) => (if s.dashpot { "DASHPOTA" } else { "SPRINGA" }, String::new(), vec![s.first, b]), + (Some(b), Some(d), false) => (if s.dashpot { "DASHPOT2" } else { "SPRING2" }, format!("{d}, {d}"), vec![s.first, b]), + (None, Some(d), _) => (if s.dashpot { "DASHPOT1" } else { "SPRING1" }, d.to_string(), vec![s.first]), + (_, None, _) => { + bad.insert(format!("{what} `{}`: a direction along no global axis (needs a local orientation, not written yet)", s.name)); + continue; + } + }; + push(&s.name, what, element_type, DiscreteProperty::Spring { dashpot: s.dashpot, dofs, value: s.value }, (label()?, nodes)); + } + unsupported(FORMAT, bad)?; + Ok(groups) +} + /// Append a real number in the shortest form that reads back to the same value: `0.1`, /// `210000.0`, `7.85E-9`, `1.0E300`. fn push_number(out: &mut String, v: f64) -> Result<()> { @@ -342,15 +444,7 @@ fn coupling_surfaces(model: &Model, r: &Resolved) -> BTreeMap { r.sets.iter().filter(|(n, s)| !s.faces.is_empty() || !r.planar_edges(n).is_empty()).map(|(n, _)| n.to_ascii_uppercase()).collect(); let mut out = BTreeMap::new(); for name in model.couplings().keys() { - let base = format!("{name}-nodes"); - let mut candidate = base.clone(); - let mut k = 2u32; - while taken.contains(&candidate.to_ascii_uppercase()) { - candidate = format!("{base}-{k}"); - k = k.saturating_add(1); - } - taken.insert(candidate.to_ascii_uppercase()); - out.insert(name.clone(), candidate); + out.insert(name.clone(), crate::unique_name(&format!("{name}-nodes"), &mut taken)); } out } @@ -581,6 +675,41 @@ impl Deck { } } + fn discrete_elements(&mut self, groups: &[DiscreteGroup]) { + if groups.is_empty() { + return; + } + self.comment("Point masses, springs and dashpots"); + for g in groups { + self.line(&format!("*ELEMENT, TYPE={}, ELSET={}", g.element_type, g.elset)); + for (label, nodes) in &g.members { + let nodes: Vec = nodes.iter().map(u32::to_string).collect(); + self.line(&format!("{label}, {}", nodes.join(", "))); + } + } + } + + fn discrete_properties(&mut self, groups: &[DiscreteGroup]) -> Result<()> { + for g in groups { + match &g.property { + DiscreteProperty::Mass(m) => { + self.line(&format!("*MASS, ELSET={}", g.elset)); + self.numbers_line("", &[*m])?; + } + DiscreteProperty::Inertia(i) => { + self.line(&format!("*ROTARY INERTIA, ELSET={}", g.elset)); + self.numbers_line("", i)?; + } + DiscreteProperty::Spring { dashpot, dofs, value } => { + self.line(&format!("{}, ELSET={}", if *dashpot { "*DASHPOT" } else { "*SPRING" }, g.elset)); + self.line(if dofs.is_empty() { "," } else { dofs }); + self.numbers_line("", &[*value])?; + } + } + } + Ok(()) + } + fn sets(&mut self, model: &Model, r: &Resolved) -> Result<()> { if r.sets.is_empty() { return Ok(()); diff --git a/crates/decks/src/lib.rs b/crates/decks/src/lib.rs index 076d236..973aaa9 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 modelcraft_model::nodal::edge_pressure; -use modelcraft_model::{Diagnostic, Load, Model, Resolved}; +use modelcraft_model::{Diagnostic, Load, Model, Resolved, ResolvedDirection}; #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum Error { @@ -158,3 +158,28 @@ pub fn edge_pressure_loads<'a>(model: &Model, r: &Resolved, active: impl IntoIte out.retain(|_, v| *v != 0.0); Ok(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 { + match d { + ResolvedDirection::Dof(k) => Some(k), + ResolvedDirection::Vector(n) | ResolvedDirection::Line(n) => { + let axis = n.iter().position(|c| c.abs() == 1.0)?; + (n.iter().filter(|c| **c != 0.0).count() == 1).then_some(axis as u8 + 1) + } + } +} + +/// `base`, or `base-2`, `base-3`, … : the first that is not in `taken` (compared without case), +/// which it is then added to. +pub(crate) fn unique_name(base: &str, taken: &mut BTreeSet) -> String { + let mut candidate = base.to_string(); + let mut k = 2u32; + while taken.contains(&candidate.to_ascii_uppercase()) { + candidate = format!("{base}-{k}"); + k = k.saturating_add(1); + } + taken.insert(candidate.to_ascii_uppercase()); + candidate +} diff --git a/crates/decks/src/lsdyna/elements.rs b/crates/decks/src/lsdyna/elements.rs index 980dbc6..855edef 100644 --- a/crates/decks/src/lsdyna/elements.rs +++ b/crates/decks/src/lsdyna/elements.rs @@ -38,7 +38,7 @@ impl Keyword { ElementKind::Hex20 => Keyword::SolidH20, ElementKind::Tet4 | ElementKind::Tet10 | ElementKind::Wedge6 | ElementKind::Hex8 => Keyword::Solid, ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => Keyword::Shell, - ElementKind::Rod2 | ElementKind::Beam2 => Keyword::Beam, + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => Keyword::Beam, } } @@ -56,6 +56,7 @@ impl Keyword { pub(crate) fn order(kind: ElementKind) -> &'static [usize] { match kind { ElementKind::Rod2 | ElementKind::Beam2 => &[0, 1], + ElementKind::Rod3 => &[0, 1, 2], ElementKind::Tri3 => &[0, 1, 2, 2], ElementKind::Tri6 => &[0, 1, 2, 2, 3, 4, 5, 5], ElementKind::Quad4 => &[0, 1, 2, 3], @@ -141,6 +142,7 @@ pub(crate) fn elform(kind: ElementKind, section: &SectionKind, f: Formulation, s }, SectionKind::Rod { .. } => match kind { ElementKind::Rod2 => plain(3), + ElementKind::Rod3 => no("LS-DYNA's trusses have two nodes; there is no counterpart for"), _ => no("a rod section cannot hold"), }, SectionKind::Beam { .. } => no("beam sections are not written to LS-DYNA decks yet (planned with beam profiles):"), diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 8e6cd6a..c011f97 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -31,6 +31,15 @@ //! (with `*LOAD_BODY_PARTS` when it acts on a set); output as `*DATABASE_BINARY_D3PLOT`, //! `*DATABASE_SPCFORC` and `*DATABASE_EXTENT_BINARY`. //! +//! Point masses are `*ELEMENT_MASS` (one mass for every direction) or `*ELEMENT_INERTIA` (the +//! inertia tensor's components, with the mass); springs and dashpots are `*ELEMENT_DISCRETE` with +//! one `*PART`, `*SECTION_DISCRETE` (torsional for a rotational degree of freedom) and +//! `*MAT_SPRING_ELASTIC` or `*MAT_DAMPER_VISCOUS` per definition, acting along the line of their +//! nodes or along a `*DEFINE_SD_ORIENTATION` vector. A spring or dashpot to the ground connects +//! to an extra node at the same place, held by `*BOUNDARY_SPC_NODE` (labels after the mesh's). +//! Element, part, section and material ids continue after the mesh's. A mass that differs by +//! direction and the three-node truss (rod3) are reported. +//! //! Several steps are independent analyses in the model (nothing carries over), which is what //! LS-DYNA's `*CASE` runs: each step is a case between `*CASE_BEGIN_n` and `*CASE_END_n`, run in //! order when LS-DYNA is started with `CASE` on its command line. A single step is written @@ -54,7 +63,7 @@ use std::collections::{BTreeMap, BTreeSet}; use modelcraft_model::{ BoundaryCondition, CouplingKind, Diagnostic, Dofs, Elastic, Element, ElementKind, Hardening, Load, Material, Model, OutputVariable, Plastic, - Procedure, Resolved, Section, SectionKind, Step, Target, + Procedure, Resolved, ResolvedDirection, Section, SectionKind, Step, Target, }; use crate::{Error, Result}; @@ -161,6 +170,7 @@ fn write_with(model: &Model, r: &Resolved) -> Result { element_cards(&mut d, &plan)?; sets(&mut d, model, r, &plan)?; couplings(&mut d, model, r, &plan)?; + discrete(&mut d, model, r, &plan)?; if let Some(parts) = &plan.gravity_parts { d.line("*SET_PART_LIST"); d.header(&["sid"]); @@ -312,6 +322,16 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result> { None } }; + if r.is_planar() { + for name in model.sections().iter().filter(|(_, s)| matches!(s.kind, SectionKind::Rod { .. })).map(|(n, _)| n) { + refused.insert(format!("section `{name}`: trusses in a 2D model (LS-DYNA's 2D elements take no truss)")); + } + } + for m in &r.masses { + if m.mass[0] != m.mass[1] || m.mass[0] != m.mass[2] { + refused.insert(format!("point mass `{}`: a mass that differs by direction (LS-DYNA's mass elements take one mass)", m.name)); + } + } let mut curves = BTreeMap::new(); for name in materials.keys() { let Some(m) = model.material(name) else { continue }; @@ -646,6 +666,122 @@ fn couplings(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan) -> Result<( Ok(()) } +/// Point masses, springs and dashpots (see the module docs). +fn discrete(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan) -> Result<()> { + if r.masses.is_empty() && r.springs.is_empty() { + return Ok(()); + } + let limit = |what: &str| Error::Limit(format!("no {what} id left for point masses, springs and dashpots")); + let mut eid = model.elements().keys().next_back().copied().unwrap_or(0); + let mut next_eid = || -> Result { + eid = eid.checked_add(1).filter(|e| *e <= modelcraft_model::MAX_LABEL).ok_or_else(|| limit("element"))?; + Ok(eid) + }; + d.comment("Point masses"); + for m in &r.masses { + let mass = m.mass[0]; + if m.inertia.iter().all(|c| *c == 0.0) { + d.line("*ELEMENT_MASS"); + d.header(&["eid", "nid", "mass"]); + d.card([Field::from(next_eid()?), Field::from(m.node), Field::from(mass)]); + } else { + let [i11, i22, i33, i12, i13, i23] = m.inertia; + d.line("*ELEMENT_INERTIA"); + d.header(&["eid", "nid", "csid"]); + d.card([Field::from(next_eid()?), Field::from(m.node)]); + d.header(&["ixx", "ixy", "ixz", "iyy", "iyz", "izz", "mass"]); + d.card([i11, i12, i13, i22, i23, i33, mass].map(Field::from)); + } + } + if r.springs.is_empty() { + return Ok(()); + } + // Ids after the mesh's parts, the rigid bodies' parts and the materials. + let kinematic = model.couplings().values().filter(|c| c.kind == CouplingKind::Kinematic).count(); + let mut pid = u32::try_from(plan.parts.len() + kinematic).map_err(|_| limit("part"))?; + let mut mid = u32::try_from(plan.materials.len()).map_err(|_| limit("material"))?; + let mut vid = 0u32; + // Ground nodes: labels after the mesh's, skipping the reference points'. + let taken: BTreeSet = model.reference_points().values().map(|p| p.label).collect(); + let mut node = model.nodes().keys().next_back().copied().unwrap_or(0); + let mut ground: Vec<(u32, [f64; 3])> = Vec::new(); + let mut elements: Vec> = Vec::new(); + let mut current: Option<&str> = None; + let mut vid_of_part = 0u32; + for s in &r.springs { + let rotational = s.direction.is_rotational(); + if current != Some(s.name.as_str()) { + current = Some(&s.name); + pid = pid.checked_add(1).ok_or_else(|| limit("part"))?; + mid = mid.checked_add(1).ok_or_else(|| limit("material"))?; + d.comment(&format!("{} `{}`", s.kind(), s.name)); + if s.dashpot { + d.line("*MAT_DAMPER_VISCOUS_TITLE"); + d.line(&s.name); + d.header(&["mid", "dc"]); + } else { + d.line("*MAT_SPRING_ELASTIC_TITLE"); + d.line(&s.name); + d.header(&["mid", "k"]); + } + d.card([Field::from(mid), Field::from(s.value)]); + d.line("*SECTION_DISCRETE"); + d.header(&["secid", "dro"]); + d.card([Field::from(pid), Field::from(i32::from(rotational))]); + d.line("*PART"); + d.line(&s.name); + d.header(&["pid", "secid", "mid"]); + d.card([Field::from(pid), Field::from(pid), Field::from(mid)]); + let axis = (!matches!(s.direction, ResolvedDirection::Line(_))).then(|| s.direction.axis()); + let v = match axis { + Some([x, y, z]) => { + vid += 1; + d.line("*DEFINE_SD_ORIENTATION"); + d.header(&["vid", "iop", "xt", "yt", "zt"]); + d.card([Field::from(vid), Field::from(0), Field::from(x), Field::from(y), Field::from(z)]); + vid + } + None => 0, + }; + vid_of_part = v; + } + let second = match s.second { + Some(b) => b, + None => { + node = node.checked_add(1).ok_or_else(|| limit("node"))?; + while taken.contains(&node) { + node = node.checked_add(1).ok_or_else(|| limit("node"))?; + } + if node > modelcraft_model::MAX_LABEL { + return Err(limit("node")); + } + ground.push((node, model.label_xyz(s.first).ok_or_else(|| broken(format!("node {}", s.first)))?)); + node + } + }; + elements.push(vec![Field::from(next_eid()?), Field::from(pid), Field::from(s.first), Field::from(second), Field::from(vid_of_part)]); + } + d.line("*ELEMENT_DISCRETE"); + d.header(&["eid", "pid", "n1", "n2", "vid"]); + for e in elements { + d.card(e); + } + if !ground.is_empty() { + d.comment("The ground ends of springs and dashpots, held in every direction"); + d.line("*NODE"); + d.header(&["nid", "x", "y", "z"]); + for (l, [x, y, z]) in &ground { + d.card([Field::from(*l), Field::from(*x), Field::from(*y), Field::from(*z)]); + } + d.line("*BOUNDARY_SPC_NODE"); + d.header(&["nid", "cid", "dofx", "dofy", "dofz", "dofrx", "dofry", "dofrz"]); + for (l, _) in &ground { + d.card([Field::from(*l), Field::Blank].into_iter().chain([1; 6].map(Field::from))); + } + } + Ok(()) +} + /// LS-DYNA's degree-of-freedom code for prescribed motion and nodal loads: 1–3 translations, /// 5–7 rotations. fn dof_code(dof: u8) -> i32 { diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs index 297f2a7..3455abd 100644 --- a/crates/decks/src/nastran/mod.rs +++ b/crates/decks/src/nastran/mod.rs @@ -35,7 +35,13 @@ //! (`PPLANE`/`PLPLANE` with `CPLSTN*`/`CQUADX`, `CTRIAX6`, `PAXI`) differ between Nastran //! dialects, so a deck could not be read the same everywhere; //! - centrifugal loads: `RFORCE` about an axis through the origin of the basic system (speed in -//! revolutions per time unit). +//! revolutions per time unit); +//! - point masses: `CONM1`, the node's 6 × 6 mass matrix (the masses per direction on the +//! translations, the inertia tensor's components on the rotations, so no product-of-inertia sign +//! convention is involved); springs and dashpots along a global degree of freedom: `CELAS2` and +//! `CDAMP2` (the second grid left blank for one to the ground). Their ids follow the rigid +//! elements'. A spring or dashpot along a skew direction and the three-node truss (rod3) have +//! no counterpart here and are reported. //! //! Field widths, continuations and number forms are in [`card`]. Model content Nastran cannot //! express is reported as [`Error::Unsupported`], every item at once, and nothing is written. @@ -96,7 +102,7 @@ pub fn order(kind: ElementKind) -> &'static [usize] { /// The connection entry's name for an element kind. pub fn element_card(kind: ElementKind) -> &'static str { match kind { - ElementKind::Rod2 => "CROD", + ElementKind::Rod2 | ElementKind::Rod3 => "CROD", ElementKind::Beam2 => "CBEAM", ElementKind::Tri3 => "CTRIA3", ElementKind::Tri6 => "CTRIA6", @@ -257,6 +263,8 @@ struct Deck<'m> { element_pid: BTreeMap, /// Rigid elements: (id, coupling name, coupling). rigid: Vec<(u32, &'m str, &'m Coupling)>, + /// Point masses, springs and dashpots: (comment, card). + discrete: Vec<(String, Card)>, /// Node → (dependent degrees of freedom, coupling that makes them dependent). dependent: BTreeMap, /// Constraint sets by id. @@ -304,6 +312,7 @@ impl<'m> Deck<'m> { properties: Vec::new(), element_pid: BTreeMap::new(), rigid: Vec::new(), + discrete: Vec::new(), dependent: BTreeMap::new(), spc_sets: Vec::new(), load_sids: BTreeMap::new(), @@ -315,6 +324,7 @@ impl<'m> Deck<'m> { self.plan_materials()?; self.plan_sections()?; self.plan_couplings()?; + self.plan_discrete()?; self.plan_constraints()?; self.plan_loads()?; for sub in &mut self.subcases { @@ -447,6 +457,9 @@ impl<'m> Deck<'m> { self.add_property(sname, Property::Beam { mid, profile, orientation }, elements.iter().copied())?; } SectionKind::Rod { area } => { + if elements.iter().any(|e| kind_of(*e) == Some(ElementKind::Rod3)) { + self.problems.insert(format!("section `{sname}`: three-node trusses (rod3; Nastran's rods have two grids)")); + } let Some(mid) = ids.isotropic else { continue }; self.add_property(sname, Property::Rod { mid, area }, elements.iter().copied())?; } @@ -520,6 +533,44 @@ impl<'m> Deck<'m> { Ok(()) } + /// Point masses, springs and dashpots, numbered after the rigid elements. + fn plan_discrete(&mut self) -> Result<()> { + let r = self.r; + let mut next = self.rigid.last().map(|e| e.0).or_else(|| self.model.elements().keys().next_back().copied()).unwrap_or(0); + let mut id = || -> Result { + next = next + .checked_add(1) + .filter(|n| *n <= MAX_LABEL) + .ok_or_else(|| Error::Limit(format!("no element id left below {} for point masses, springs and dashpots", MAX_LABEL as u64 + 1)))?; + Ok(next) + }; + for m in &r.masses { + let [mx, my, mz] = m.mass; + let [i11, i22, i33, i12, i13, i23] = m.inertia; + // M11, M21, M22, M31, M32, M33, then the rotations' rows: M41…M44, M51…M55, M61…M66. + let card = Card::new("CONM1") + .int(id()?) + .int(m.node) + .blank() + .reals([mx, 0.0, my, 0.0, 0.0, mz, 0.0, 0.0, 0.0, i11, 0.0, 0.0, 0.0, i12, i22, 0.0, 0.0, 0.0, i13, i23, i33]); + self.discrete.push((format!("point mass {}", m.name), card)); + } + for s in &r.springs { + let what = s.kind(); + let Some(dof) = crate::axis_dof(s.direction) else { + self.problems.insert(format!("{what} `{}`: a direction along no global axis (CELAS and CDAMP act along grid components)", s.name)); + continue; + }; + let card = Card::new(if s.dashpot { "CDAMP2" } else { "CELAS2" }).int(id()?).real(s.value).int(s.first).int(dof); + let card = match s.second { + Some(b) => card.int(b).int(dof), + None => card, + }; + self.discrete.push((format!("{what} {}", s.name), card)); + } + Ok(()) + } + /// One constraint set per distinct step constraint, checked against coupled nodes. fn plan_constraints(&mut self) -> Result<()> { let (model, r) = (self.model, self.r); @@ -780,6 +831,13 @@ impl<'m> Deck<'m> { }; write_named(out, &format!("coupling {name}"), *eid, card)?; } + if !self.discrete.is_empty() { + comment(out, "point masses, springs and dashpots"); + } + for (what, card) in &self.discrete { + comment(out, what); + card.write(out)?; + } Ok(()) } diff --git a/crates/decks/tests/discrete.rs b/crates/decks/tests/discrete.rs new file mode 100644 index 0000000..1c2aa3e --- /dev/null +++ b/crates/decks/tests/discrete.rs @@ -0,0 +1,166 @@ +//! Trusses (rod3), point masses, rotary inertia, springs and dashpots in every deck format: +//! what each writes, and what each reports as unsupported. + +use modelcraft_decks::{Error, abaqus, lsdyna, nastran}; +use modelcraft_model::{ + BoundaryCondition, Direction, Elastic, ElementKind, Formulation, Mass, Material, Model, PointMass, Procedure, ReferencePoint, Section, + SectionKind, Set, Spring, Step, Target, UnitSystem, +}; + +/// A rod2 (or rod3) bar from node 1 to node 2 and a rod2 on to node 4, a reference point beyond +/// it; a motor mass at the tip, a hub with inertia on the point, a spring to the ground across +/// the tip, a spring along the line from the tip to the hub, a dashpot to the ground. +fn model(rod3: bool) -> Model { + let mut m = Model::default(); + m.set_units(UnitSystem::MmTonneS); + for (l, x) in [(1, 0.0), (2, 1000.0), (3, 500.0), (4, 2000.0)] { + m.set_node(l, [x, 0.0, 0.0]).unwrap(); + } + if rod3 { + m.set_element(1, ElementKind::Rod3, vec![1, 2, 3]).unwrap(); + } else { + m.set_element(1, ElementKind::Rod2, vec![1, 2]).unwrap(); + } + m.set_element(2, ElementKind::Rod2, vec![2, 4]).unwrap(); + m.add_set("bars", Set { elements: vec![1, 2], ..Set::default() }).unwrap(); + m.add_set("base", Set { nodes: vec![1], ..Set::default() }).unwrap(); + m.add_set("tip", Set { nodes: vec![4], ..Set::default() }).unwrap(); + m.add_material( + "Steel", + Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() }, + ) + .unwrap(); + m.add_section( + "bars", + Section { set: "bars".into(), material: "Steel".into(), kind: SectionKind::Rod { area: 20.0 }, formulation: Formulation::default() }, + ) + .unwrap(); + m.add_reference_point("Hub", ReferencePoint { xyz: [2500.0, 0.0, 0.0], label: 0 }).unwrap(); + let tip = Target::Set("tip".into()); + let hub = Target::ReferencePoint("Hub".into()); + m.add_mass("motor", PointMass { target: tip.clone(), mass: Some(Mass::Isotropic(0.002)), inertia: None }).unwrap(); + m.add_mass("hub", PointMass { target: hub.clone(), mass: Some(Mass::Isotropic(0.01)), inertia: Some([1.0, 2.0, 3.0, 0.1, 0.0, 0.0]) }).unwrap(); + m.add_spring("mount", Spring::Spring { first: tip.clone(), second: None, direction: Direction::Dof(2), stiffness: 50.0 }).unwrap(); + m.add_spring("link", Spring::Spring { first: tip.clone(), second: Some(hub), direction: Direction::Line, stiffness: 1000.0 }).unwrap(); + m.add_spring("damper", Spring::Dashpot { first: tip, second: None, direction: Direction::Vector([0.0, 0.0, -2.0]), coefficient: 0.5 }).unwrap(); + m.add_boundary_condition("fix", BoundaryCondition::Fixed { target: Target::Set("base".into()) }).unwrap(); + let mut step = Step::new("Modes", Procedure::Frequency { modes: 3, min_hz: None, max_hz: None }); + step.boundary_conditions.insert("fix".into()); + m.add_step(step, None).unwrap(); + m +} + +fn unsupported(e: std::result::Result) -> String { + match e { + Err(Error::Unsupported { what, .. }) => what, + other => panic!("expected unsupported, got {other:?}"), + } +} + +#[test] +fn abaqus_writes_every_discrete_element_and_rod3() { + let deck = abaqus::write(&model(true)).unwrap(); + for expected in [ + "*ELEMENT, TYPE=T3D3\n1, 1, 2, 3\n", + "*ELEMENT, TYPE=MASS, ELSET=hub-mass\n3, 99999999\n", + "*ELEMENT, TYPE=MASS, ELSET=motor-mass\n4, 4\n", + "*ELEMENT, TYPE=ROTARYI, ELSET=hub-inertia\n5, 99999999\n", + "*ELEMENT, TYPE=DASHPOT1, ELSET=damper-dashpot\n6, 4\n", + "*ELEMENT, TYPE=SPRINGA, ELSET=link-spring\n7, 4, 99999999\n", + "*ELEMENT, TYPE=SPRING1, ELSET=mount-spring\n8, 4\n", + "*MASS, ELSET=hub-mass\n0.01\n", + "*ROTARY INERTIA, ELSET=hub-inertia\n1.0, 2.0, 3.0, 0.1, 0.0, 0.0\n", + "*DASHPOT, ELSET=damper-dashpot\n3\n0.5\n", + "*SPRING, ELSET=link-spring\n,\n1000.0\n", + "*SPRING, ELSET=mount-spring\n2\n50.0\n", + ] { + assert!(deck.contains(expected), "{expected}\n{deck}"); + } + // An element set name that a model set already has is made unique. + let mut m = model(false); + m.add_set("motor-mass", Set { nodes: vec![4], ..Set::default() }).unwrap(); + assert!(abaqus::write(&m).unwrap().contains("ELSET=motor-mass-2\n")); + // Skew directions and masses per direction need orientations. + let mut m = model(false); + m.replace_spring( + "mount", + Spring::Spring { first: Target::Set("tip".into()), second: None, direction: Direction::Vector([0.0, 1.0, 1.0]), stiffness: 1.0 }, + ) + .unwrap(); + m.replace_mass("motor", PointMass { target: Target::Set("tip".into()), mass: Some(Mass::Directional([1.0, 2.0, 3.0])), inertia: None }).unwrap(); + let what = unsupported(abaqus::write(&m)); + assert!(what.contains("spring `mount`") && what.contains("point mass `motor`"), "{what}"); +} + +#[test] +fn nastran_writes_conm1_celas2_cdamp2() { + let deck = nastran::write(&model(false)).unwrap(); + eprintln!("{deck}"); + for expected in ["CONM1*", "CELAS2*", "CDAMP2*", "$ point mass hub", "$ spring link", "$ dashpot damper"] { + assert!(deck.contains(expected), "{expected}\n{deck}"); + } + // The hub's 6 × 6 mass matrix: mass on the translations, the inertia tensor on the rotations. + let conm1: Vec<&str> = deck.lines().skip_while(|l| !l.starts_with("CONM1*")).take(6).collect(); + let fields: Vec = conm1 + .iter() + .flat_map(|l| l.get(8..).unwrap_or("").as_bytes().chunks(16).map(|c| std::str::from_utf8(c).unwrap().trim().to_string())) + .filter(|t| !t.is_empty()) + .filter_map(|t| t.parse().ok()) + .collect(); + assert_eq!(fields.first().copied(), Some(3.0), "element id after the mesh's: {conm1:?}"); + assert!( + fields.ends_with(&[0.01, 0.0, 0.01, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.1, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0]), + "{fields:?}" + ); + // The line spring runs along X: component 1 at both grids; the dashpot along −Z: component 3. + let link = deck.lines().skip_while(|l| !l.contains("spring link")).nth(1).unwrap().to_string(); + assert!(link.starts_with("CELAS2*") && link.contains("1000.") && link.ends_with("4 1"), "{link}"); + let rest = deck.lines().skip_while(|l| !l.contains("spring link")).nth(2).unwrap().to_string(); + assert!(rest.contains("99999999") && rest.trim_end().ends_with('1'), "{rest}"); + let what = unsupported(nastran::write(&model(true))); + assert!(what.contains("rod3"), "{what}"); + let mut m = model(false); + m.replace_spring( + "link", + Spring::Spring { first: Target::Set("tip".into()), second: None, direction: Direction::Vector([1.0, 1.0, 0.0]), stiffness: 1.0 }, + ) + .unwrap(); + assert!(unsupported(nastran::write(&m)).contains("spring `link`")); +} + +#[test] +fn lsdyna_writes_masses_inertia_and_discrete_springs() { + let deck = lsdyna::write(&model(false)).unwrap().text; + eprintln!("{deck}"); + for expected in [ + "*ELEMENT_MASS", + "*ELEMENT_INERTIA", + "*MAT_SPRING_ELASTIC_TITLE\nlink\n", + "*MAT_DAMPER_VISCOUS_TITLE\ndamper\n", + "*SECTION_DISCRETE", + "*DEFINE_SD_ORIENTATION", + "*ELEMENT_DISCRETE", + "*BOUNDARY_SPC_NODE", + ] { + assert!(deck.contains(expected), "{expected}\n{deck}"); + } + // Two grounded links (damper, mount) get two ground nodes after the mesh's (5, 6); the + // line spring has no orientation vector (vid 0, along its nodes). + let discrete: Vec> = deck + .lines() + .skip_while(|l| *l != "*ELEMENT_DISCRETE") + .skip(1) + .take_while(|l| !l.starts_with('*')) + .filter(|l| !l.starts_with('$')) + .map(|l| l.split_whitespace().map(|f| f.parse().unwrap()).collect()) + .collect(); + assert_eq!(discrete.len(), 3); + assert_eq!(discrete[0][2..4], [4, 5], "damper to its ground node"); + assert_eq!(discrete[1][2..4], [4, 99_999_999], "link to the hub"); + assert_eq!(discrete[1].get(4).copied().unwrap_or(0), 0, "along the line"); + assert_eq!(discrete[2][2..4], [4, 6], "mount to its ground node"); + assert!(unsupported(lsdyna::write(&model(true))).contains("rod3")); + let mut m = model(false); + m.replace_mass("motor", PointMass { target: Target::Set("tip".into()), mass: Some(Mass::Directional([1.0, 2.0, 3.0])), inertia: None }).unwrap(); + assert!(unsupported(lsdyna::write(&m)).contains("point mass `motor`")); +} diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index caac4a0..abef15d 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -306,8 +306,22 @@ pub static CATALOGUE: &[Capability] = &[ soon("interaction.contact.find", "Find Contacts", Interactions, "Contact", "M7", "Detect contact pairs automatically"), soon("interaction.rigid_wall", "Rigid Wall", Interactions, "Contact", "M10", "Rigid walls for explicit dynamics"), soon("interaction.bolt", "Bolt Load", Interactions, "Connectors", "M7", "Bolt pretension"), - soon("interaction.spring", "Spring", Interactions, "Connectors", "M7", "Springs and dashpots"), - soon("interaction.mass", "Point Mass", Interactions, "Connectors", "M7", "Point masses and inertia"), + ready_more( + "model.spring.create", + "Spring", + Interactions, + "Connectors", + "Spring or dashpot to the ground or between two points, along a direction or their line", + &["model.spring.edit", "model.spring.delete"], + ), + ready_more( + "model.mass.create", + "Point Mass", + Interactions, + "Connectors", + "Point mass and rotary inertia at nodes or a reference point", + &["model.mass.edit", "model.mass.delete"], + ), soon("interaction.connector", "Connector", Interactions, "Connectors", "M7", "Connector elements between points"), soon("interaction.fsi", "FSI Interface", Interactions, "Coupled", "M13", "The fluid–structure interface"), // ---- Solve ---- diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs index c95217a..a5c913b 100644 --- a/crates/engine/src/cmd/define.rs +++ b/crates/engine/src/cmd/define.rs @@ -1,13 +1,14 @@ //! Model definition commands: create, edit and delete sets, materials, sections, reference -//! points, couplings, boundary conditions, loads, fluid domains and steps, and the step activity -//! table. +//! points, couplings, point masses, springs and dashpots, boundary conditions, loads, fluid +//! domains and steps, and the step activity table. //! //! Each definition's parameters are its JSON form in the model (the same form a `.modelcraft` //! file stores), plus its `name`. `create` fails when the name is taken, `edit` replaces an //! existing definition, `delete` refuses while another definition refers to it. use modelcraft_model::{ - BoundaryCondition, Coupling, FaceCondition, FluidBoundary, FluidDomain, Load, Material, Model, ReferencePoint, Section, Set, Step, + BoundaryCondition, Coupling, FaceCondition, FluidBoundary, FluidDomain, Load, Material, Model, PointMass, ReferencePoint, Section, Set, Spring, + Step, }; use serde::de::DeserializeOwned; use serde_json::{Value, json}; @@ -176,6 +177,34 @@ definition_commands!( COUPLING_DELETE ); +definition_commands!( + "point mass", + PointMass, + "mass", + "Point Mass", + r#"{"name": "motor", "target": {"set": "…"} | {"reference_point": "…"}, "mass": 2.5 | [mx, my, mz] (per global direction), "inertia": [I11, I22, I33, I12, I13, I23] (tensor about the node in global axes; needs rotations: a reference point)} — at every node of the target, each with the full value"#, + add_mass, + replace_mass, + remove_mass, + MASS_CREATE, + MASS_EDIT, + MASS_DELETE +); + +definition_commands!( + "spring", + Spring, + "spring", + "Spring or Dashpot", + r#"{"name": "mount", "type": "spring", "stiffness": 1000 | "type": "dashpot", "coefficient": 5, "first": {"set": "…"} | {"reference_point": "…"}, "second": (optional; without it each node of `first` connects to the ground; with it each end is one node), "direction": {"dof": 1..6} (4–6: torsional) | {"vector": [x, y, z]} | "line" (along the two ends, two-point only)} — dashpots act in dynamic steps; static and frequency steps ignore them"#, + add_spring, + replace_spring, + remove_spring, + SPRING_CREATE, + SPRING_EDIT, + SPRING_DELETE +); + definition_commands!( "boundary condition", BoundaryCondition, diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index f23d4eb..bc01e75 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -54,6 +54,12 @@ pub static ALL: &[CommandSpec] = &[ define::COUPLING_CREATE, define::COUPLING_EDIT, define::COUPLING_DELETE, + define::MASS_CREATE, + define::MASS_EDIT, + define::MASS_DELETE, + define::SPRING_CREATE, + define::SPRING_EDIT, + define::SPRING_DELETE, define::BOUNDARY_CONDITION_CREATE, define::BOUNDARY_CONDITION_EDIT, define::BOUNDARY_CONDITION_DELETE, diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index 54ce8da..f27423d 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -950,4 +950,32 @@ pub(crate) mod tests { // Defining steps edited the model, so the impact result is stale. assert_eq!(s.solved_steps().count(), 0); } + + /// A mass on a spring and a dashpot, defined, solved and written through commands; editing + /// and deleting them goes through the same checks as every definition. + #[test] + fn spring_mass_through_commands() { + let mut s = Session::default(); + let script = json!([ + {"command": "mesh.node.add", "params": {"label": 1, "xyz": [0, 0, 0]}}, + {"command": "model.set.create", "params": {"name": "n", "nodes": [1]}}, + {"command": "model.mass.create", "params": {"name": "m", "target": {"set": "n"}, "mass": 4}}, + {"command": "model.spring.create", "params": {"name": "k", "type": "spring", "first": {"set": "n"}, "direction": {"dof": 1}, "stiffness": 400}}, + {"command": "model.spring.create", "params": {"name": "c", "type": "dashpot", "first": {"set": "n"}, "direction": {"vector": [1, 0, 0]}, "coefficient": 1}}, + {"command": "model.boundary_condition.create", "params": {"name": "yz", "type": "displacement", "target": {"set": "n"}, "values": [null, 0, 0, null, null, null]}}, + {"command": "model.step.create", "params": {"name": "Modes", "procedure": {"type": "frequency", "modes": 1}, "boundary_conditions": ["yz"]}}, + ]); + s.run_script(&script).unwrap(); + let r = s.execute("solve.frequency", &json!({})).unwrap(); + let f = r["modes"][0]["frequency"].as_f64().unwrap(); + assert!((f - 10.0 / (2.0 * std::f64::consts::PI)).abs() < 1e-12, "{r}"); + assert!(r["warnings"].as_array().unwrap().iter().any(|w| w.as_str().unwrap().contains("dashpots `c`")), "{r}"); + assert!(s.execute("model.mass.edit", &json!({"name": "m", "target": {"set": "n"}, "mass": -1})).is_err()); + assert!(s.execute("model.set.delete", &json!({"name": "n"})).is_err(), "the set is in use"); + s.execute("model.mass.edit", &json!({"name": "m", "target": {"set": "n"}, "mass": 1})).unwrap(); + let r = s.execute("solve.frequency", &json!({})).unwrap(); + assert!((r["modes"][0]["frequency"].as_f64().unwrap() - 20.0 / (2.0 * std::f64::consts::PI)).abs() < 1e-12); + s.execute("model.spring.delete", &json!({"name": "c"})).unwrap(); + assert!(s.model.springs().len() == 1 && s.model.mass("m").is_some()); + } } diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index c25fc74..9773773 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -404,6 +404,13 @@ mod tests { json!({"name": "Bend", "procedure": {"type": "static"}, "loads": "push"}), json!({"step": "Bend", "loads": ["nope"], "active": "yes"}), json!({"step": "Bend", "boundary_conditions": [], "loads": []}), + json!({"name": "A", "target": {"set": "end"}, "mass": [1e308, -1, f64::MAX], "inertia": [1, 1, 1, 1e300, 0, 0]}), + json!({"name": "A", "target": {"reference_point": "RP"}, "inertia": [1, 2]}), + json!({"name": "A", "target": {"set": "end"}, "mass": "heavy"}), + json!({"name": "A", "type": "spring", "first": {"set": "end"}, "second": {"set": "end"}, "direction": "line", "stiffness": 1e308}), + json!({"name": "A", "type": "spring", "first": {"set": "end"}, "direction": {"dof": 255}, "stiffness": 1}), + json!({"name": "A", "type": "dashpot", "first": {"set": "end"}, "direction": {"vector": [1e-320, 0, 0]}, "coefficient": -1}), + json!({"name": "A", "type": "dashpot", "first": {"reference_point": "RP"}, "direction": {"dof": 4}, "coefficient": 1e300}), json!({"system": "furlongs"}), // Fluid domains, boundaries and steps. json!({"name": "A", "material": "Steel", "min": [0, 0, 0], "max": [1e308, 1, 1], "spacing": 1e-300}), diff --git a/crates/explicit/src/discrete.rs b/crates/explicit/src/discrete.rs new file mode 100644 index 0000000..51034b1 --- /dev/null +++ b/crates/explicit/src/discrete.rs @@ -0,0 +1,316 @@ +//! Trusses, point masses, springs and dashpots in explicit dynamics. +//! +//! - **Trusses** (rod2, rod3), elastic, with lumped (HRZ) mass. With nonlinear geometry a rod2 +//! is corotational: its strain is (l − L)/L from its current length l and its force acts along +//! its current axis (Belytschko et al. 2014, §2.4 and §4.5); a rod3 then is not supported yet. +//! Without it the internal force is K u with the kernel's stiffness ([`modelcraft_fem::truss`]). +//! The stress is the uniaxial tensor along the axis (the mean of the integration points'), and +//! the critical increment comes from the kernel's frequency bound (exact for rod2: L/c). +//! - **Point masses** add to the lumped mass of their node, per direction. +//! - **Springs** add the force k (cᵀ Δu) c and **dashpots** c_d (cᵀ Δv) c, with Δu, Δv the +//! relative motion of the ends along the unit direction c (a fixed direction; one along the line +//! of the ends follows that line as it turns under nonlinear geometry, with the stretch +//! l − L₀); dashpots use the velocity of the +//! previous half increment, v^{n−½}, as the bulk viscosity does (Belytschko, Liu, Moran, +//! Elkhodary 2014, §6.2). Spring work counts as internal energy, dashpot work as viscous +//! (dissipated) energy. +//! - **Stable increment** of springs and dashpots, per degree of freedom: with +//! ω² ≤ Σⱼ |Kᵢⱼ| / √(mᵢ mⱼ) (Gershgorin's bound on M^{-1/2} K M^{-1/2}, exact for one spring on +//! one mass) and the same row sum D of the damping, the central-difference scheme with lagged +//! damping is stable for Δt ≤ (2/ω)(√(1 + ξ²) − ξ) with ξ = D / (2ω) (Belytschko et al. 2014, +//! §6.3), which is 2/D for a dashpot alone. Rotary inertia and rotational springs need +//! rotational degrees of freedom, which this solver does not have yet. + +use modelcraft_fem::DofKind; +use modelcraft_fem::dense::Matrix; +use modelcraft_fem::truss::{Truss, uniaxial}; +use modelcraft_model::{ElementKind, Resolved, ResolvedDirection}; + +use crate::system::NODE_DOFS; +use crate::{Error, Result}; + +/// A truss element. +#[derive(Debug, Clone)] +pub struct Rod { + pub label: u32, + /// Position in [`crate::system::System::elements`]. + pub index: usize, + pub nodes: Vec, + pub material: usize, + pub young: f64, + /// Stiffness on the element's translations, node by node. + pub stiffness: Matrix, + /// Per integration point: the axial strain operator and the tangent. + pub strain: Vec<(Vec, [f64; 3])>, + /// The lumped mass per node for unit density. + pub nodal_volume: Vec, + /// The stable length: the critical increment is `length / wave_speed`. + pub length: f64, + /// E A and the rest length (a rod2's corotational force). + pub axial: f64, + pub rest_length: f64, +} + +/// What a truss is made of: Young's modulus, density and cross-sectional area. +#[derive(Debug, Clone, Copy)] +pub struct RodSection { + pub young: f64, + pub density: f64, + pub area: f64, +} + +impl Rod { + /// `material` is the index in [`crate::system::System::materials`]. + pub fn new( + label: u32, + index: usize, + kind: ElementKind, + nodes: Vec, + xyz: &[[f64; 3]], + material: usize, + section: RodSection, + ) -> Result { + let RodSection { young, density, area } = section; + let at = |e| Error::element(label, e); + let t = Truss::new(kind, xyz).map_err(at)?; + let omega2 = t.frequency_bound(young, density).map_err(at)?; + let wave_speed = (young / density).sqrt(); + Ok(Rod { + label, + index, + nodes, + material, + young, + stiffness: t.stiffness(young, area).map_err(at)?, + strain: t.strain_operator().map_err(at)?, + nodal_volume: t.lumped_mass(1.0, area).map_err(at)?, + length: 2.0 * wave_speed / omega2.sqrt(), + axial: young * area, + rest_length: t.length().map_err(at)?, + }) + } + + /// The element's degrees of freedom, node by node. + pub fn dofs(&self) -> impl Iterator + '_ { + self.nodes.iter().flat_map(|n| modelcraft_fem::truss::DOFS.node.iter().filter_map(DofKind::direction).map(|d| dof(*n, d))) + } + + /// Add its internal force into `force` and return the element's mean stress and strain + /// tensors: corotational for a rod2 with nonlinear geometry (`xyz` the rest positions of the + /// system's nodes), K u otherwise. + pub fn forces(&self, xyz: &[[f64; 3]], u: &[f64], nonlinear: bool, force: &mut [f64]) -> ([f64; 6], [f64; 6]) { + if nonlinear && let [a, b] = self.nodes[..] { + let at = |n: usize| -> [f64; 3] { + let x = xyz.get(n).copied().unwrap_or_default(); + [0, 1, 2].map(|k| x[k] + u.get(dof(n, k)).copied().unwrap_or(0.0)) + }; + let (pa, pb) = (at(a), at(b)); + let d = [pb[0] - pa[0], pb[1] - pa[1], pb[2] - pa[2]]; + let l = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt(); + let t = if l > 0.0 { d.map(|c| c / l) } else { [0.0; 3] }; + let e = (l - self.rest_length) / self.rest_length; + let n = self.axial * e; + for k in 0..3 { + if let Some(x) = force.get_mut(dof(a, k)) { + *x -= n * t[k]; + } + if let Some(x) = force.get_mut(dof(b, k)) { + *x += n * t[k]; + } + } + return (uniaxial(self.young * e, t, 1.0), uniaxial(e, t, 2.0)); + } + // At most three nodes of three translations: the element vectors live on the stack. + let n = self.stiffness.n.min(MAX_DOFS); + let mut ue = [0.0; MAX_DOFS]; + for (x, d) in ue.iter_mut().zip(self.dofs()) { + *x = u.get(d).copied().unwrap_or(0.0); + } + for (i, d) in self.dofs().enumerate().take(n) { + let f: f64 = (0..n).map(|j| self.stiffness.get(i, j) * ue[j]).sum(); + if let Some(x) = force.get_mut(d) { + *x += f; + } + } + let (mut stress, mut strain) = ([0.0; 6], [0.0; 6]); + let share = 1.0 / self.strain.len().max(1) as f64; + for (row, t) in &self.strain { + let e: f64 = row.iter().zip(&ue).map(|(a, b)| a * b).sum(); + for (k, (s, x)) in uniaxial(self.young * e, *t, 1.0).iter().zip(uniaxial(e, *t, 2.0)).enumerate() { + stress[k] += share * s; + strain[k] += share * x; + } + } + (stress, strain) + } +} + +/// Degrees of freedom of the largest truss (rod3: three nodes of three translations). +const MAX_DOFS: usize = 9; + +/// The global degree of freedom of translation `d` (0–2) of node `n` in the system's uniform, +/// node-major layout ([`NODE_DOFS`] per node, checked when the system is built). +pub fn dof(n: usize, d: usize) -> usize { + NODE_DOFS * n + d +} + +/// The node and the direction (0–2) of a degree of freedom: the inverse of [`dof`]. +pub fn node_direction(dof: usize) -> (usize, usize) { + (dof / NODE_DOFS, dof % NODE_DOFS) +} + +/// A spring or dashpot on degrees of freedom: the elongation is Σ c u(second) − Σ c u(first). +#[derive(Debug, Clone)] +pub struct Link { + pub name: String, + pub dashpot: bool, + /// (degree of freedom, coefficient) with the first end's coefficients negated. + pub terms: Vec<(usize, f64)>, + pub value: f64, + /// For one along the line of its ends: the two nodes and their rest distance. + pub line: Option<(usize, usize, f64)>, +} + +impl Link { + fn stretch(&self, x: &[f64]) -> f64 { + self.terms.iter().map(|(d, c)| c * x.get(*d).copied().unwrap_or(0.0)).sum() + } + + /// Add its internal force (spring on `u`, dashpot on `v`; `xyz` the rest positions of the + /// system's nodes, for a line that turns under nonlinear geometry). + pub fn forces(&self, xyz: &[[f64; 3]], u: &[f64], v: &[f64], nonlinear: bool, stiffness: &mut [f64], viscous: &mut [f64]) { + if nonlinear && let Some((a, b, rest)) = self.line { + let at = |n: usize, x: &[f64]| -> [f64; 3] { [0, 1, 2].map(|k| x.get(dof(n, k)).copied().unwrap_or(0.0)) }; + let pos = |n: usize| { + let (p, d) = (xyz.get(n).copied().unwrap_or_default(), at(n, u)); + [p[0] + d[0], p[1] + d[1], p[2] + d[2]] + }; + let (pa, pb) = (pos(a), pos(b)); + let d = [pb[0] - pa[0], pb[1] - pa[1], pb[2] - pa[2]]; + let l = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt(); + let t = if l > 0.0 { d.map(|c| c / l) } else { [0.0; 3] }; + let (va, vb) = (at(a, v), at(b, v)); + let s = if self.dashpot { self.value * (0..3).map(|k| t[k] * (vb[k] - va[k])).sum::() } else { self.value * (l - rest) }; + let out = if self.dashpot { viscous } else { stiffness }; + for k in 0..3 { + if let Some(x) = out.get_mut(dof(a, k)) { + *x -= s * t[k]; + } + if let Some(x) = out.get_mut(dof(b, k)) { + *x += s * t[k]; + } + } + return; + } + let (x, out) = if self.dashpot { (v, viscous) } else { (u, stiffness) }; + let s = self.value * self.stretch(x); + for (d, c) in &self.terms { + if let Some(f) = out.get_mut(*d) { + *f += s * c; + } + } + } +} + +/// Springs, dashpots and point masses of a model on the explicit system's nodes. +#[derive(Debug, Clone, Default)] +pub struct Discrete { + pub links: Vec, + /// Point mass per degree of freedom. + pub mass: Vec, +} + +impl Discrete { + /// `node` maps a label to its node index. + /// `node` maps a label to its node index, `xyz` gives the nodes' rest positions. + pub fn new(resolved: &Resolved, dofs: usize, xyz: &[[f64; 3]], node: impl Fn(u32) -> Option) -> Result { + let off = |what: &str, name: &str, label: u32| { + Error::NotSupported(format!( + "{what} `{name}` at node {label}: explicit dynamics takes discrete elements on mesh nodes (not reference points) so far" + )) + }; + let mut out = Discrete { links: Vec::new(), mass: vec![0.0; dofs] }; + for m in &resolved.masses { + if m.inertia.iter().any(|c| *c != 0.0) { + return Err(Error::NotSupported(format!("rotary inertia (point mass `{}`) in explicit dynamics: it has no rotations yet", m.name))); + } + let n = node(m.node).ok_or_else(|| off("point mass", &m.name, m.node))?; + for d in 0..3 { + if let Some(x) = out.mass.get_mut(dof(n, d)) { + *x += m.mass[d]; + } + } + } + for s in &resolved.springs { + let what = s.kind(); + if s.direction.is_rotational() { + return Err(Error::NotSupported(format!("rotational {what}s (`{}`) in explicit dynamics: it has no rotations yet", s.name))); + } + let c = s.direction.coefficients(); + let mut terms = Vec::new(); + let mut ends = Vec::new(); + for (label, sign) in std::iter::once((s.first, -1.0)).chain(s.second.map(|b| (b, 1.0))) { + let n = node(label).ok_or_else(|| off(what, &s.name, label))?; + ends.push(n); + terms.extend((0..3).filter(|d| c[*d] != 0.0).map(|d| (dof(n, d), sign * c[d]))); + } + let line = match (s.direction, ends.as_slice()) { + (ResolvedDirection::Line(_), [a, b]) => { + let (pa, pb) = (xyz.get(*a).copied().unwrap_or_default(), xyz.get(*b).copied().unwrap_or_default()); + Some((*a, *b, ((pb[0] - pa[0]).powi(2) + (pb[1] - pa[1]).powi(2) + (pb[2] - pa[2]).powi(2)).sqrt())) + } + _ => None, + }; + out.links.push(Link { name: s.name.clone(), dashpot: s.dashpot, terms, value: s.value, line }); + } + Ok(out) + } + + /// The smallest stable increment of the springs and dashpots on the lumped mass `mass` (per + /// degree of freedom), and the name of a link at the degree of freedom that sets it. A link + /// along its turning line (`nonlinear`) is bounded for every direction it may take. + pub fn critical_increment(&self, mass: &[f64], nonlinear: bool) -> Option<(f64, String)> { + let mut stiff = vec![0.0; mass.len()]; + let mut damp = vec![0.0; mass.len()]; + let mut owner: Vec> = vec![None; mass.len()]; + for (k, l) in self.links.iter().enumerate() { + let rows = if l.dashpot { &mut damp } else { &mut stiff }; + let turning: Vec<(usize, f64)>; + let terms = match l.line { + Some((a, b, _)) if nonlinear => { + turning = (0..3).map(|k| (dof(a, k), -1.0)).chain((0..3).map(|k| (dof(b, k), 1.0))).collect(); + &turning + } + _ => &l.terms, + }; + for (i, ci) in terms { + let mi = mass.get(*i).copied().unwrap_or(0.0); + let sum: f64 = terms.iter().map(|(j, cj)| (l.value * ci * cj).abs() / (mi * mass.get(*j).copied().unwrap_or(0.0)).sqrt()).sum(); + if let Some(r) = rows.get_mut(*i) { + *r += sum; + } + if let Some(o) = owner.get_mut(*i) { + o.get_or_insert(k); + } + } + } + let mut best: Option<(f64, usize)> = None; + for (i, (w2, d)) in stiff.iter().zip(&damp).enumerate() { + let dt = match (*w2 > 0.0, *d > 0.0) { + (false, false) => continue, + (false, true) => 2.0 / d, + (true, _) => { + let w = w2.sqrt(); + let xi = d / (2.0 * w); + (2.0 / w) * ((1.0 + xi * xi).sqrt() - xi) + } + }; + if best.is_none_or(|(b, _)| dt < b) { + best = Some((dt, i)); + } + } + let (dt, i) = best?; + let name = owner.get(i).copied().flatten().and_then(|k| self.links.get(k)).map(|l| l.name.clone()).unwrap_or_default(); + Some((dt, name)) + } +} diff --git a/crates/explicit/src/elements.rs b/crates/explicit/src/elements.rs index 211eb2e..c2e3b2f 100644 --- a/crates/explicit/src/elements.rs +++ b/crates/explicit/src/elements.rs @@ -477,6 +477,22 @@ pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Wor critical = critical.min(element.finish(ws, &es.point.stress, &inc.v, m.bulk)); } } + for r in &sys.rods { + let (sigma, strain) = r.forces(&sys.xyz, m.u, sys.nonlinear, &mut ws.stress_force); + if let Some(es) = state.elements.get_mut(r.index) { + es.point.stress = sigma; + es.strain = strain; + } + if let Some(mat) = sys.materials.get(r.material) { + critical = critical.min(r.length / mat.wave_speed); + } + } + for l in &sys.discrete.links { + l.forces(&sys.xyz, m.u, m.v, sys.nonlinear, &mut ws.stress_force, &mut ws.viscous_force); + } + if let Some((dt, _)) = &sys.discrete_increment { + critical = critical.min(*dt); + } Ok(critical) } diff --git a/crates/explicit/src/integrate.rs b/crates/explicit/src/integrate.rs index 9c78e50..8f9a9b0 100644 --- a/crates/explicit/src/integrate.rs +++ b/crates/explicit/src/integrate.rs @@ -133,13 +133,31 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res .iter() .map(|h| (h.label, h.length, h.material)) .chain(sys.tets.iter().map(|t| (t.label, t.length, t.material))) + .chain(sys.rods.iter().map(|r| (r.label, r.length, r.material))) .map(|(label, length, m)| ControllingElement { label, increment: length / sys.materials.get(m).map_or(f64::INFINITY, |m| m.wave_speed) }) .collect(); critical.sort_by(|a, b| a.increment.total_cmp(&b.increment).then(a.label.cmp(&b.label))); critical.truncate(CONTROLLING); - let smallest = critical - .first() - .map(|c| c.increment) + let n = sys.dofs(); + let mass = sys.dof_mass(); + if let Some(d) = mass.iter().position(|m| !(*m > 0.0 && m.is_finite())) { + let label = sys.labels.get(crate::discrete::node_direction(d).0).copied().unwrap_or(0); + return Err(Error::Invalid(format!( + "node {label} has no mass in direction {}; explicit dynamics needs a mass at every node (give the nodes only springs or dashpots use a point mass)", + crate::discrete::node_direction(d).1 + 1 + ))); + } + // Springs and dashpots have their own stable increment on the lumped mass. + let elements_smallest = critical.first().map(|c| c.increment); + if let Some((l, name)) = &sys.discrete_increment + && elements_smallest.is_none_or(|e| *l < e) + { + warnings.push(format!("spring or dashpot `{name}` limits the time increment to {l:.3e} (before the scale factor)")); + } + let smallest = elements_smallest + .into_iter() + .chain(sys.discrete_increment.as_ref().map(|l| l.0)) + .reduce(f64::min) .filter(|d| *d > 0.0 && d.is_finite()) .ok_or_else(|| Error::Invalid("no element has a finite, positive critical time increment (check the material constants)".into()))?; let estimate = end / (scale * smallest); @@ -170,12 +188,8 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res ))); } - let n = sys.dofs(); - let mass: Vec = (0..n).map(|d| sys.mass.get(d / NODE_DOFS).copied().unwrap_or(0.0)).collect(); - if mass.iter().any(|m| !(*m > 0.0 && m.is_finite())) { - return Err(Error::Internal("a node has no positive mass".into())); - } - let total_mass: f64 = sys.mass.iter().sum(); + // Point masses count with their mean over the three directions. + let total_mass: f64 = sys.mass.iter().sum::() + sys.discrete.mass.iter().sum::() / 3.0; let mut ws = Workspace::new(&sys); let mut previous = Workspace::new(&sys); let mut state = State::new(&sys); @@ -211,17 +225,17 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res history.time.push(t); history.increment.push(dt); history.energies.push(e); - let sum = |x: &[f64], weight: &dyn Fn(usize) -> f64| { + let sum = |x: &[f64], weight: &dyn Fn(usize, usize) -> f64| { let mut s = [0.0; 3]; for (node, c) in x.as_chunks::().0.iter().enumerate() { for k in 0..3 { - s[k] += weight(node) * c[k]; + s[k] += weight(node, k) * c[k]; } } s }; - history.momentum.push(sum(v, &|node| sys.mass.get(node).copied().unwrap_or(0.0))); - history.constraint_force.push(sum(r, &|_| 1.0)); + history.momentum.push(sum(v, &|node, k| mass.get(crate::discrete::dof(node, k)).copied().unwrap_or(0.0))); + history.constraint_force.push(sum(r, &|_, _| 1.0)); let pick = |x: &[f64]| sys.history_nodes.iter().map(|i| [0, 1, 2].map(|k| x.get(NODE_DOFS * i + k).copied().unwrap_or(0.0))).collect(); history.displacement.push(pick(u)); history.velocity.push(pick(v)); diff --git a/crates/explicit/src/lib.rs b/crates/explicit/src/lib.rs index c6dbfb0..75e9f8d 100644 --- a/crates/explicit/src/lib.rs +++ b/crates/explicit/src/lib.rs @@ -32,6 +32,7 @@ // Small fixed-size vector kernels read as the textbook formulas with index loops. #![allow(clippy::needless_range_loop)] +pub mod discrete; pub mod elements; pub mod integrate; pub mod result; diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs index e326e4f..b2b5428 100644 --- a/crates/explicit/src/system.rs +++ b/crates/explicit/src/system.rs @@ -1,5 +1,6 @@ -//! From the resolved model to the explicit system: elements with their materials and stable -//! lengths, lumped masses, prescribed velocities, rigid walls, initial velocities and the load +//! From the resolved model to the explicit system: elements (solids and trusses) with their +//! materials and stable lengths, lumped masses with the point masses, springs and dashpots +//! ([`crate::discrete`]), prescribed velocities, rigid walls, initial velocities and the load //! vector of an explicit dynamic step. use modelcraft_fem::material::Elasticity; @@ -10,6 +11,7 @@ use modelcraft_model::{ Elastic, ElementKind, ExplicitDynamic, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step, StressRate, Target, }; +use crate::discrete::{Discrete, Rod}; use crate::elements::{HexScales, add_at, frequency_bound, hex_scales}; use crate::{Error, Result}; @@ -97,8 +99,14 @@ pub struct System { pub materials: Vec, pub hexes: Vec, pub tets: Vec, + pub rods: Vec, /// Row-sum lumped mass per node. pub mass: Vec, + /// Point masses, springs and dashpots. + pub discrete: Discrete, + /// The springs' and dashpots' stable increment on the lumped mass, and a link that sets it, + /// when there are any. + pub discrete_increment: Option<(f64, String)>, /// Prescribed velocity per degree of freedom (`None` = free). pub prescribed: Vec>, pub walls: Vec, @@ -162,6 +170,7 @@ fn check_formulation(kind: ElementKind, section_name: &str, f: &modelcraft_model ElementKind::Hex8 => Ok(()), ElementKind::Tet4 if f.integration == Integration::Reduced => refuse("reduced integration of tet4 elements"), ElementKind::Tet4 => Ok(()), + ElementKind::Rod2 | ElementKind::Rod3 => Ok(()), other => refuse(&format!("{} elements in explicit dynamics (it has one-point hex8 and nodal-pressure tet4 so far)", other.name())), } } @@ -173,14 +182,18 @@ impl System { } // Elements by label, with their section's material. let mut materials: Vec = Vec::new(); - let mut chosen: Vec<(u32, ElementKind, Vec, usize)> = Vec::new(); + let mut chosen: Vec<(u32, ElementKind, Vec, usize, Option)> = Vec::new(); for (name, section) in model.sections() { - if !matches!(section.kind, SectionKind::Solid {}) { - return Err(Error::NotSupported(format!( - "{} sections (section `{name}`); explicit dynamics has solid elements so far", - section.kind.name() - ))); - } + let area = match section.kind { + SectionKind::Solid {} => None, + SectionKind::Rod { area } => Some(area), + _ => { + return Err(Error::NotSupported(format!( + "{} sections (section `{name}`); explicit dynamics has solid and truss elements so far", + section.kind.name() + ))); + } + }; let mi = match materials.iter().position(|m| m.name == section.material) { Some(i) => i, None => { @@ -191,6 +204,13 @@ impl System { let set = resolved.sets.get(§ion.set).ok_or_else(|| invalid(format!("set `{}` is not resolved", section.set)))?; for &label in &set.elements { let el = model.element(label).ok_or_else(|| invalid(format!("element {label} is missing")))?; + if area.is_some() != matches!(el.kind, ElementKind::Rod2 | ElementKind::Rod3) { + return Err(invalid(format!( + "section `{name}`: a {} section cannot hold {} element {label}", + section.kind.name(), + el.kind.name() + ))); + } check_formulation(el.kind, name, §ion.formulation)?; if el.kind == ElementKind::Tet4 && materials.get(mi).is_some_and(|m| m.j2.is_none()) { return Err(Error::NotSupported(format!( @@ -198,11 +218,13 @@ impl System { section.material ))); } - chosen.push((label, el.kind, el.nodes.clone(), mi)); + chosen.push((label, el.kind, el.nodes.clone(), mi, area)); } } chosen.sort_by_key(|e| e.0); let mut labels: Vec = chosen.iter().flat_map(|e| e.2.iter().copied()).collect(); + // Mesh nodes only discrete elements use (reference points are refused by Discrete::new). + labels.extend(resolved.discrete_labels().filter(|l| model.node(*l).is_some())); labels.sort_unstable(); labels.dedup(); let layout = DofLayout::uniform(labels.iter().copied(), DofSet::TRANSLATIONS); @@ -213,7 +235,7 @@ impl System { return Err(Error::Internal("explicit dynamics needs three translations per node".into())); } if labels.is_empty() { - return Err(invalid("the model has no elements to analyse")); + return Err(invalid("the model has no elements or point masses to analyse")); } let xyz: Vec<[f64; 3]> = labels.iter().map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| invalid(format!("node {l} is missing")))).collect::>()?; @@ -224,7 +246,10 @@ impl System { materials, hexes: Vec::new(), tets: Vec::new(), + rods: Vec::new(), mass: vec![0.0; n], + discrete: Discrete::default(), + discrete_increment: None, prescribed: vec![None; dofs], walls: Vec::new(), initial_velocity: vec![0.0; dofs], @@ -238,7 +263,7 @@ impl System { labels, xyz, }; - for (index, (label, kind, nodes, mi)) in chosen.into_iter().enumerate() { + for (index, (label, kind, nodes, mi, area)) in chosen.into_iter().enumerate() { let idx: Vec = nodes.iter().map(|l| sys.node(*l).ok_or_else(|| Error::Internal(format!("node {l}")))).collect::>()?; let mat = sys.materials.get(mi).ok_or_else(|| Error::Internal("material index".into()))?; let (rho, c) = (mat.density, mat.wave_speed); @@ -264,9 +289,34 @@ impl System { sys.lump(&nodes, &mass); sys.tets.push(Tet { label, index, nodes, material: mi, kinematics: k, mass, omega2, length: 2.0 * c / omega2.sqrt() }); } + ElementKind::Rod2 | ElementKind::Rod3 => { + let area = area.ok_or_else(|| Error::Internal("truss area".into()))?; + let source = model.material(&mat.name).ok_or_else(|| Error::Internal("truss material".into()))?; + if source.plastic.is_some() || source.damage.is_some() { + return Err(Error::NotSupported(format!( + "plastic or damaging materials on trusses (material `{}`, element {label}): explicit trusses are elastic so far", + mat.name + ))); + } + if kind == ElementKind::Rod3 && settings.nonlinear_geometry { + return Err(Error::NotSupported(format!( + "three-node trusses with nonlinear geometry (element {label}); set the step's nonlinear_geometry to false" + ))); + } + let young = + modelcraft_fem::truss::young_modulus(source.elastic.as_ref().ok_or_else(|| Error::Internal("truss material".into()))?) + .map_err(|e| invalid(format!("material `{}` of element {label}: {e}", mat.name)))?; + let x: Vec<[f64; 3]> = idx.iter().map(|i| sys.xyz.get(*i).copied().unwrap_or_default()).collect(); + let rod = Rod::new(label, index, kind, idx.clone(), &x, mi, crate::discrete::RodSection { young, density: rho, area })?; + let mass: Vec = rod.nodal_volume.iter().map(|v| rho * v).collect(); + sys.lump(&idx, &mass); + sys.rods.push(rod); + } other => return Err(Error::Internal(format!("{} element kind", other.name()))), } } + sys.discrete = Discrete::new(resolved, sys.dofs(), &sys.xyz, |l| sys.node(l))?; + sys.discrete_increment = sys.discrete.critical_increment(&sys.dof_mass(), sys.nonlinear); sys.prescribe(model, resolved, step, settings.end_time)?; sys.load(model, resolved, step)?; sys.conditions(resolved, settings)?; @@ -290,6 +340,13 @@ impl System { self.layout.dofs() } + /// The lumped mass per degree of freedom: the elements' and the point masses'. + pub fn dof_mass(&self) -> Vec { + (0..self.dofs()) + .map(|d| self.mass.get(crate::discrete::node_direction(d).0).copied().unwrap_or(0.0) + self.discrete.mass.get(d).copied().unwrap_or(0.0)) + .collect() + } + /// 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<()> { @@ -363,6 +420,12 @@ impl System { add_at(&mut f, *node, &acceleration.map(|a| rho * w * a)); } }); + // Point masses take gravity on the whole model. + if set.is_none() { + for (d, (x, m)) in f.iter_mut().zip(&self.discrete.mass).enumerate() { + *x += m * acceleration.get(crate::discrete::node_direction(d).1).copied().unwrap_or(0.0); + } + } } Load::BodyForce { set, vector } => { let chosen = resolved.sets.get(set).map(|r| r.elements.as_slice()).unwrap_or_default(); @@ -399,6 +462,9 @@ impl System { for t in self.tets.iter().filter(|t| wanted(t.label)) { visit(&t.nodes, &t.kinematics.nodal_volume, density(t.material)); } + for r in self.rods.iter().filter(|r| wanted(r.label)) { + visit(&r.nodes, &r.nodal_volume, density(r.material)); + } } /// The analysed nodes of a resolved set. diff --git a/crates/explicit/tests/one_d.rs b/crates/explicit/tests/one_d.rs new file mode 100644 index 0000000..01ae8b8 --- /dev/null +++ b/crates/explicit/tests/one_d.rs @@ -0,0 +1,307 @@ +//! Trusses, point masses, springs and dashpots in explicit dynamics against closed-form +//! solutions. Run with `--nocapture` to see the errors. +//! +//! - A mass on a spring started with a velocity: x = (v₀/ω) sin ωt. +//! - The same with a dashpot: the amplitude decays as e^{−ξωt}, the period is 2π/ω_d with +//! ω_d = ω√(1 − ξ²); a mass on a dashpot alone slows as v₀ e^{−ct/m}. +//! - Two masses on a two-point spring along its line: ω² = k(1/m₁ + 1/m₂), momentum conserved. +//! - A truss bar driven at one end: a wave front at c t with stress −ρ c v behind it. + +use std::f64::consts::PI; + +use modelcraft_explicit::{ExplicitResult, Limits, solve_explicit}; +use modelcraft_fem::truss::axial; +use modelcraft_model::*; + +fn model(step: ExplicitDynamic) -> Model { + let mut m = Model::default(); + m.add_step(Step::new("Dynamic", Procedure::ExplicitDynamic(ExplicitDynamic::new(step.end_time))), None).unwrap(); + // The step's settings name sets; they are put in once the sets exist (see `solve`). + PENDING.with(|p| *p.borrow_mut() = Some(step)); + m +} + +thread_local! { + static PENDING: std::cell::RefCell> = const { std::cell::RefCell::new(None) }; +} + +/// Put the step's settings in, keeping what is active in it, and solve. +fn solve(m: &mut Model) -> ExplicitResult { + if let Some(settings) = PENDING.with(|p| p.borrow_mut().take()) { + let mut step = m.step("Dynamic").unwrap().clone(); + step.procedure = Procedure::ExplicitDynamic(settings); + m.replace_step(step).unwrap(); + } + solve_explicit(m, None, &Limits::default()).unwrap() +} + +fn set(m: &mut Model, name: &str, nodes: &[u32]) -> Target { + m.add_set(name, Set { nodes: nodes.to_vec(), ..Set::default() }).unwrap(); + Target::Set(name.into()) +} + +fn hold(m: &mut Model, name: &str, target: Target, values: [Option; 6]) { + m.add_boundary_condition(name, BoundaryCondition::Displacement { target, values }).unwrap(); + m.set_boundary_condition_active("Dynamic", name, true).unwrap(); +} + +fn step(end: f64, scale: f64, moving: &str, velocity: [f64; 3]) -> ExplicitDynamic { + let mut s = ExplicitDynamic::new(end); + s.dt_scale = scale; + s.history_nodes = Some(moving.into()); + s.history_interval = Some(end / 20_000.0); + s.initial_velocities = vec![InitialVelocity { set: moving.into(), velocity }]; + s.bulk_viscosity = BulkViscosity { linear: 0.0, quadratic: 0.0 }; + s +} + +/// The history of one direction of the first history node. +fn trace(r: &ExplicitResult, k: usize) -> Vec<(f64, f64)> { + r.history.time.iter().zip(&r.history.displacement).map(|(t, d)| (*t, d[0][k])).collect() +} + +/// Local maxima of a sampled signal (parabolic interpolation through the three samples). +fn peaks(x: &[(f64, f64)]) -> Vec<(f64, f64)> { + x.windows(3) + .filter(|w| w[1].1 > w[0].1 && w[1].1 >= w[2].1) + .map(|w| { + let (h, a, b, c) = (w[1].0 - w[0].0, w[0].1, w[1].1, w[2].1); + let d = 0.5 * (a - c) / (a - 2.0 * b + c); + (w[1].0 + d * h, b - 0.25 * (a - c) * d) + }) + .collect() +} + +/// A mass on a spring to the ground (and optionally a dashpot), started at velocity v₀. +fn oscillator(k: f64, m: f64, damping: Option, periods: f64, scale: f64) -> ExplicitResult { + let omega = (k / m).sqrt(); + let v0 = 0.3; + let mut model = model(step(periods * 2.0 * PI / omega, scale, "n", [v0, 0.0, 0.0])); + model.set_node(1, [0.0, 0.0, 0.0]).unwrap(); + let n = set(&mut model, "n", &[1]); + model.add_spring("k", Spring::Spring { first: n.clone(), second: None, direction: Direction::Dof(1), stiffness: k }).unwrap(); + if let Some(c) = damping { + model.add_spring("c", Spring::Dashpot { first: n.clone(), second: None, direction: Direction::Dof(1), coefficient: c }).unwrap(); + } + model.add_mass("m", PointMass { target: n.clone(), mass: Some(Mass::Isotropic(m)), inertia: None }).unwrap(); + hold(&mut model, "yz", n, [None, Some(0.0), Some(0.0), None, None, None]); + solve(&mut model) +} + +#[test] +fn spring_mass_oscillation() { + let (k, m) = (200.0f64, 0.5); + let r = oscillator(k, m, None, 3.0, 0.02); + let omega = (k / m).sqrt(); + let amplitude = 0.3 / omega; + let worst = trace(&r, 0).iter().map(|(t, x)| (x - amplitude * (omega * t).sin()).abs()).fold(0.0, f64::max); + let p = peaks(&trace(&r, 0)); + let period = (p[2].0 - p[0].0) / 2.0; + eprintln!( + "spring–mass: worst error {:.2e} of the amplitude, period {period:.6} vs {:.6}, Δt {:.3e}, energy error {:.1e}", + worst / amplitude, + 2.0 * PI / omega, + r.initial_increment, + r.energy.relative_error(r.initial_kinetic) + ); + assert!(worst < 2e-3 * amplitude, "second-order phase error"); + assert!((period - 2.0 * PI / omega).abs() < 1e-4 * period); + assert!(r.energy.relative_error(r.initial_kinetic) < 1e-3); + assert!(r.warnings.iter().any(|w| w.contains("limits the time increment")), "{:?}", r.warnings); + // The increment is the scale times the spring's 2/ω. + assert!((r.initial_increment - 0.02 * 2.0 / omega).abs() < 1e-12 * r.initial_increment); +} + +#[test] +fn damped_decay_rate() { + let (k, m, xi) = (200.0f64, 0.5, 0.05); + let omega = (k / m).sqrt(); + let c = 2.0 * xi * m * omega; + let r = oscillator(k, m, Some(c), 6.0, 0.01); + let p = peaks(&trace(&r, 0)); + let wd = omega * (1.0 - xi * xi).sqrt(); + // Decay rate from the first and last peaks, period from their spacing. + let (first, last) = (p[0], p[p.len() - 1]); + let rate = (first.1 / last.1).ln() / (last.0 - first.0); + let period = (last.0 - first.0) / (p.len() - 1) as f64; + let rate_err = (rate - xi * omega).abs() / (xi * omega); + let period_err = (period - 2.0 * PI / wd).abs() / (2.0 * PI / wd); + eprintln!("damped oscillator: decay rate {rate:.6} vs ξω {:.6} ({rate_err:.2e}), period error {period_err:.2e}", xi * omega); + assert!(rate_err < 0.01, "{rate} vs {}", xi * omega); + assert!(period_err < 1e-3); + // The dashpot dissipates what the oscillator loses: the balance holds. + assert!(r.energy.viscous > 0.9 * r.initial_kinetic && r.energy.relative_error(r.initial_kinetic) < 0.01, "{:?}", r.energy); + + // A dashpot alone: v = v₀ e^{−ct/m}. + let mut model = model(step(3.0 * m / c, 0.001, "n", [0.3, 0.0, 0.0])); + model.set_node(1, [0.0; 3]).unwrap(); + let n = set(&mut model, "n", &[1]); + model.add_spring("c", Spring::Dashpot { first: n.clone(), second: None, direction: Direction::Vector([2.0, 0.0, 0.0]), coefficient: c }).unwrap(); + model.add_mass("m", PointMass { target: n.clone(), mass: Some(Mass::Isotropic(m)), inertia: None }).unwrap(); + hold(&mut model, "yz", n, [None, Some(0.0), Some(0.0), None, None, None]); + let r = solve(&mut model); + let (t, v) = (r.history.time.last().unwrap(), r.history.velocity.last().unwrap()[0][0]); + let exact = 0.3 * (-c * t / m).exp(); + eprintln!("dashpot alone: v {v:.6e} vs {exact:.6e} at t = {t:.4}"); + assert!((v - exact).abs() < 0.01 * exact); + // Travel: x(∞) = v₀ m / c. + let x = r.history.displacement.last().unwrap()[0][0]; + let travel = 0.3 * m / c * (1.0 - (-c * t / m).exp()); + assert!((x - travel).abs() < 0.01 * travel, "{x} vs {travel}"); +} + +#[test] +fn two_masses_on_a_line_spring() { + let (k, m1, m2) = (50.0f64, 1.0, 3.0); + let omega = (k * (1.0 / m1 + 1.0 / m2)).sqrt(); + let dir = [0.6, 0.0, 0.8]; + let mut model = model(step(4.0 * 2.0 * PI / omega, 0.02, "a", [0.0; 3])); + model.set_node(1, [0.0; 3]).unwrap(); + model.set_node(2, dir.map(|c| 2.0 * c)).unwrap(); + let a = set(&mut model, "a", &[1]); + let b = set(&mut model, "b", &[2]); + model.add_spring("k", Spring::Spring { first: a.clone(), second: Some(b.clone()), direction: Direction::Line, stiffness: k }).unwrap(); + model.add_mass("m1", PointMass { target: a, mass: Some(Mass::Isotropic(m1)), inertia: None }).unwrap(); + model.add_mass("m2", PointMass { target: b, mass: Some(Mass::Isotropic(m2)), inertia: None }).unwrap(); + // Node 1 starts towards node 2, node 2 at rest: the centre of mass drifts, the spring rings. + PENDING.with(|p| p.borrow_mut().as_mut().unwrap().initial_velocities = vec![InitialVelocity { set: "a".into(), velocity: dir.map(|c| 0.5 * c) }]); + let r = solve(&mut model); + // Momentum is conserved exactly (internal forces are equal and opposite). + let p0 = r.history.momentum[0]; + let drift = r.history.momentum.iter().map(|p| (0..3).map(|k| (p[k] - p0[k]).abs()).fold(0.0, f64::max)).fold(0.0, f64::max); + // The relative stretch along the line oscillates at ω. + let x1 = trace(&r, 0); + let along: Vec<(f64, f64)> = r + .history + .time + .iter() + .zip(&r.history.displacement) + .map(|(t, d)| (*t, d[0][0] * dir[0] + d[0][2] * dir[2] - 0.5 * m1 / (m1 + m2) * t)) + .collect(); + let p = peaks(&along); + let period = (p[3].0 - p[0].0) / 3.0; + eprintln!("two masses: period {period:.6} vs {:.6}, momentum drift {drift:.1e}", 2.0 * PI / omega); + assert!((period - 2.0 * PI / omega).abs() < 1e-4 * period); + assert!(drift < 1e-14); + assert!(x1.iter().all(|(_, x)| x.is_finite())); +} + +/// A rod2 bar along x, held across, its end pushed at velocity v: a wave front at c t and +/// σ = −ρ c v behind it, carried at the bar's own critical increment L/c. +#[test] +fn velocity_wave_in_a_truss_bar() { + let (young, rho, area, length, cells) = (1.0e4f64, 1.0, 0.01, 10.0, 200usize); + let c = (young / rho).sqrt(); + let t = 0.6 * length / c; + let v = 0.01; + for kind in [ElementKind::Rod2, ElementKind::Rod3] { + let mut s = ExplicitDynamic::new(t); + s.output_interval = Some(t / 2.0); + // A rod3 is small-displacement only; the rod2 runs corotational (the default). + s.nonlinear_geometry = kind == ElementKind::Rod2; + let mut m = model(s); + let h = length / cells as f64; + for i in 0..=cells { + m.set_node(1 + i as u32, [i as f64 * h, 0.0, 0.0]).unwrap(); + } + for i in 0..cells { + let mut nodes = vec![1 + i as u32, 2 + i as u32]; + if kind == ElementKind::Rod3 { + m.set_node(1000 + i as u32, [(i as f64 + 0.5) * h, 0.0, 0.0]).unwrap(); + nodes.push(1000 + i as u32); + } + m.set_element(1 + i as u32, kind, nodes).unwrap(); + } + m.add_set("bar", Set { elements: (1..=cells as u32).collect(), ..Set::default() }).unwrap(); + m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young, poisson: 0.0 }), density: Some(rho), ..Material::default() }) + .unwrap(); + m.add_section( + "bar", + Section { set: "bar".into(), material: "Mat".into(), kind: SectionKind::Rod { area }, formulation: Formulation::default() }, + ) + .unwrap(); + let all: Vec = m.nodes().keys().copied().collect(); + let line = set(&mut m, "line", &all); + hold(&mut m, "line", line, [None, Some(0.0), Some(0.0), None, None, None]); + let start = set(&mut m, "start", &[1]); + hold(&mut m, "push", start, [Some(v * t), None, None, None, None, None]); + let r = solve(&mut m); + let f = r.frames.last().unwrap(); + // σ behind the front (x < 0.8 c t) and the front where it falls to half. + let stress: Vec<(f64, f64)> = (0..cells).map(|i| ((i as f64 + 0.5) * h, axial(&f.stress[i], [1.0, 0.0, 0.0], false))).collect(); + let behind: Vec = stress.iter().filter(|(x, _)| *x > 0.1 * c * t && *x < 0.8 * c * t).map(|s| s.1).collect(); + let sigma = behind.iter().sum::() / behind.len() as f64; + let front = + stress.windows(2).find(|w| w[0].1 / (-rho * c * v) >= 0.5 && w[1].1 / (-rho * c * v) < 0.5).map(|w| (w[0].0 + w[1].0) / 2.0).unwrap(); + let expected_dt = if kind == ElementKind::Rod2 { 0.9 * h / c } else { r.initial_increment }; + eprintln!( + "{kind:?} bar wave: σ {sigma:.6} vs −ρcv {:.6}, front {front:.4} vs c t {:.4}, Δt {:.3e}, energy error {:.1e}", + -rho * c * v, + c * t, + r.initial_increment, + r.energy.relative_error(r.initial_kinetic) + ); + assert!((sigma + rho * c * v).abs() < 0.01 * rho * c * v); + assert!((front - c * t).abs() < 3.0 * h); + assert!((r.initial_increment - expected_dt).abs() < 1e-12 * expected_dt); + assert!(r.energy.relative_error(r.initial_kinetic) < 0.01, "{:?}", r.energy); + } +} + +#[test] +fn what_explicit_cannot_take_is_named() { + // A node with a spring and no mass. + let mut m = model(ExplicitDynamic::new(1.0)); + m.set_node(1, [0.0; 3]).unwrap(); + let n = set(&mut m, "n", &[1]); + m.add_spring("k", Spring::Spring { first: n.clone(), second: None, direction: Direction::Dof(1), stiffness: 1.0 }).unwrap(); + let e = solve_explicit(&m, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("node 1 has no mass"), "{e}"); + // Rotary inertia, and a rotational spring. + m.add_reference_point("RP", ReferencePoint { xyz: [0.0; 3], label: 0 }).unwrap(); + m.add_mass("j", PointMass { target: Target::ReferencePoint("RP".into()), mass: None, inertia: Some([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) }).unwrap(); + let e = solve_explicit(&m, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("rotary inertia"), "{e}"); +} + +/// A mass whirling on a stiff rod2 tethered to a fixed point (nonlinear geometry): the rod turns +/// with it and keeps its length, carrying the centripetal force m v²/L. A small-displacement rod +/// would let the mass fly off along the tangent. +#[test] +fn whirling_mass_on_a_corotational_rod() { + let (l, mass, v, axial) = (1.0f64, 1.0, 1.0, 1.0e4); + let t = std::f64::consts::FRAC_PI_2 * l / v; + let mut s = step(t, 0.9, "tip", [0.0, v, 0.0]); + s.bulk_viscosity = BulkViscosity { linear: 0.0, quadratic: 0.0 }; + let mut m = model(s); + m.set_node(1, [0.0; 3]).unwrap(); + m.set_node(2, [l, 0.0, 0.0]).unwrap(); + m.set_element(1, ElementKind::Rod2, vec![1, 2]).unwrap(); + m.add_set("rod", Set { elements: vec![1], ..Set::default() }).unwrap(); + // A light rod (E A = 1e4, ρ A = 1e-4) with the mass at its end. + m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: axial, poisson: 0.0 }), density: Some(1e-4), ..Material::default() }) + .unwrap(); + m.add_section( + "rod", + Section { set: "rod".into(), material: "Mat".into(), kind: SectionKind::Rod { area: 1.0 }, formulation: Formulation::default() }, + ) + .unwrap(); + let pivot = set(&mut m, "pivot", &[1]); + m.add_boundary_condition("pivot", BoundaryCondition::Fixed { target: pivot }).unwrap(); + m.set_boundary_condition_active("Dynamic", "pivot", true).unwrap(); + let tip = set(&mut m, "tip", &[2]); + m.add_mass("m", PointMass { target: tip.clone(), mass: Some(Mass::Isotropic(mass)), inertia: None }).unwrap(); + hold(&mut m, "plane", tip, [None, None, Some(0.0), None, None, None]); + let r = solve(&mut m); + let radius = r.history.displacement.iter().map(|d| ((l + d[0][0]).powi(2) + d[0][1].powi(2)).sqrt()).fold(0.0f64, |w, x| w.max((x - l).abs())); + let end = r.history.displacement.last().unwrap()[0]; + eprintln!( + "whirling mass: largest change of radius {radius:.3e} (centripetal stretch m v²/EA = {:.1e}), end at ({:.4}, {:.4})", + mass * v * v / axial, + l + end[0], + end[1] + ); + assert!(radius < 5.0 * mass * v * v / axial, "the rod keeps its length to its elastic stretch"); + // A quarter turn: the mass ends near (0, L). + assert!((l + end[0]).abs() < 1e-2 && (end[1] - l).abs() < 1e-2, "{end:?}"); +} diff --git a/crates/fem/src/basis.rs b/crates/fem/src/basis.rs index 5ee7c47..ff1d6eb 100644 --- a/crates/fem/src/basis.rs +++ b/crates/fem/src/basis.rs @@ -199,6 +199,7 @@ impl Basis { pub fn of(kind: ElementKind) -> Basis { match kind { ElementKind::Rod2 | ElementKind::Beam2 => Basis::Line2, + ElementKind::Rod3 => Basis::Line3, ElementKind::Tri3 => Basis::Tri3, ElementKind::Tri6 => Basis::Tri6, ElementKind::Quad4 => Basis::Quad4, diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs index f9decff..42f8a1d 100644 --- a/crates/fem/src/lib.rs +++ b/crates/fem/src/lib.rs @@ -24,6 +24,8 @@ //! Green–Naghdi); //! - [`plasticity`]: J2 radial return with isotropic, kinematic and Johnson–Cook hardening, //! Johnson–Cook damage; +//! - [`truss`]: two- and three-node trusses (straight or curved), their stiffness, mass, body +//! force and axial stress; //! - [`dense`]: the small dense linear algebra inside an element. //! //! Conventions: element matrices are ordered node by node, each node's degrees of freedom in @@ -46,6 +48,7 @@ pub mod plasticity; pub mod quadrature; pub mod rotation; pub mod shape; +pub mod truss; pub mod uniform; pub use basis::Basis; diff --git a/crates/fem/src/truss.rs b/crates/fem/src/truss.rs new file mode 100644 index 0000000..0c11b60 --- /dev/null +++ b/crates/fem/src/truss.rs @@ -0,0 +1,404 @@ +//! Truss (axial) elements: the two-node rod2 and the three-node rod3, straight or curved, in any +//! orientation (a plane truss is a space truss whose nodes lie in a plane). +//! +//! **Formulation.** The isoparametric bar of R. D. Cook, D. S. Malkus, M. E. Plesha, R. J. Witt, +//! *Concepts and Applications of Finite Element Analysis*, 4th ed., Wiley 2002, §3.6 and §6.2 +//! (and T. J. R. Hughes, *The Finite Element Method*, Dover 2000, §1.10), written along the +//! member's arc in three dimensions. On ξ ∈ [−1, 1] the shape functions are +//! N = ((1 − ξ)/2, (1 + ξ)/2) for rod2 and N = (ξ(ξ − 1)/2, ξ(ξ + 1)/2, 1 − ξ²) for rod3 (the +//! ends, then the middle node). The position is x(ξ) = Σ Nᵢ xᵢ, the arc length element +//! ds = J dξ with J = |x′(ξ)|, and the unit tangent t = x′/J. The axial strain is the stretch +//! along the tangent, ε = t · ∂u/∂s = Σᵢ (Nᵢ′/J) t · uᵢ, so with Bᵢ = (Nᵢ′/J) tᵀ: +//! +//! - stiffness K = ∫ E A Bᵀ B ds, +//! - consistent mass ρ A ∫ Nᵢ Nⱼ ds per direction, lumped mass by HRZ diagonal scaling (E. Hinton, +//! T. Rock, O. C. Zienkiewicz, *A note on mass lumping and related processes in the finite +//! element method*, Earthq. Eng. Struct. Dyn. 4(3):245–249, 1976), which gives ρAL/2 at each +//! end of rod2 and ρAL·(1/6, 1/6, 2/3) on a straight rod3, +//! - body force A ∫ Nᵢ b ds for a force b per volume, +//! - stress σ = E ε along the tangent: in global axes the uniaxial tensor σ t tᵀ (strain ε t tᵀ, +//! with engineering shear strains), so stresses of trusses and solids read the same way and +//! the axial force is N = A tᵀ σ t. +//! +//! The rules are Gauss–Legendre with two points for rod2 and three for rod3: exact for the +//! stiffness and mass of a straight member, and the usual full integration of a curved one. +//! Material: the Young's modulus of an isotropic material; the truss knows no other constant. + +use modelcraft_model::{Elastic, ElementKind}; + +use crate::basis::Basis; +use crate::dense::Matrix; +use crate::dof::{DofSet, ElementDofs}; +use crate::element::PointStress; +use crate::quadrature::gauss_legendre; +use crate::{Error, Result}; + +/// What a truss's nodes carry: the three translations (a plane truss is a space truss whose +/// out-of-plane motion is held; trusses of planar models will declare [`DofSet::PLANE`]). +pub const DOFS: ElementDofs = ElementDofs::new(DofSet::TRANSLATIONS); + +/// The Young's modulus a truss uses, or why there is none. +pub fn young_modulus(elastic: &Elastic) -> Result { + match *elastic { + Elastic::Isotropic { young, .. } if young > 0.0 && young.is_finite() => Ok(young), + Elastic::Isotropic { young, .. } => Err(Error::Invalid(format!("Young's modulus must be positive, got {young}"))), + Elastic::EngineeringConstants { .. } => { + Err(Error::NotSupported("orthotropic materials on trusses (a truss takes one Young's modulus: use an isotropic material)".into())) + } + } +} + +/// Shape functions and their ξ-derivatives ([`Basis::Line2`], [`Basis::Line3`]). +fn shape(n: usize, xi: f64) -> ([f64; 3], [f64; 3]) { + let s = if n == 2 { Basis::Line2 } else { Basis::Line3 }.eval([xi, 0.0, 0.0]); + let (mut v, mut d) = ([0.0; 3], [0.0; 3]); + for i in 0..n.min(3) { + v[i] = s.n[i]; + d[i] = s.dn[i][0]; + } + (v, d) +} + +/// The geometry at one point: shape functions, their arc-length derivatives Nᵢ′/J, the unit +/// tangent and the arc length per unit ξ. +#[derive(Debug, Clone)] +struct Point { + n: [f64; 3], + ds: [f64; 3], + tangent: [f64; 3], + jacobian: f64, + weight: f64, + xyz: [f64; 3], +} + +/// A truss element on its node coordinates (model node order), with the geometry at its +/// integration points worked out once. +#[derive(Debug, Clone)] +pub struct Truss<'a> { + pub kind: ElementKind, + xyz: &'a [[f64; 3]], + rule: Vec, +} + +impl<'a> Truss<'a> { + pub fn new(kind: ElementKind, xyz: &'a [[f64; 3]]) -> Result> { + if !matches!(kind, ElementKind::Rod2 | ElementKind::Rod3) { + return Err(Error::Invalid(format!("{} is not a truss element", kind.name()))); + } + if xyz.len() != kind.node_count() { + return Err(Error::Invalid(format!("{} needs {} nodes, got {}", kind.name(), kind.node_count(), xyz.len()))); + } + // The tangent must not vanish or turn back along the member (a middle node off its span). + let scale = chord(xyz).max(f64::MIN_POSITIVE); + let reference = at(xyz, 0.0, 0.0)?.tangent; + for xi in [-1.0, -0.5, 0.5, 1.0] { + let p = at(xyz, xi, 0.0)?; + if !(p.jacobian > 1e-9 * scale) || p.tangent.iter().zip(&reference).map(|(a, b)| a * b).sum::() <= 0.0 { + return Err(Error::Degenerate("the member folds back on itself or has no length".into())); + } + } + let rule = gauss_legendre(kind.node_count()).iter().map(|(xi, w)| at(xyz, *xi, *w)).collect::>()?; + Ok(Truss { kind, xyz, rule }) + } + + fn nodes(&self) -> usize { + self.kind.node_count() + } + + /// The integration points: two for rod2, three for rod3. + fn points(&self) -> Result<&[Point]> { + Ok(&self.rule) + } + + /// The natural coordinates of the integration points, in order. + pub fn abscissae(&self) -> Vec { + gauss_legendre(self.nodes()).iter().map(|(xi, _)| *xi).collect() + } + + /// The arc length. + pub fn length(&self) -> Result { + Ok(self.points()?.iter().map(|p| p.weight * p.jacobian).sum()) + } + + /// The stiffness (3n × 3n, translations node by node) for Young's modulus `young` and + /// cross-sectional area `area`. + pub fn stiffness(&self, young: f64, area: f64) -> Result { + let n = self.nodes(); + let mut k = Matrix::zeros(3 * n); + for p in self.points()? { + let b: Vec = (0..3 * n).map(|a| p.ds[a / 3] * p.tangent[a % 3]).collect(); + let factor = young * area * p.weight * p.jacobian; + for (i, bi) in b.iter().enumerate() { + for (j, bj) in b.iter().enumerate() { + k.add(i, j, factor * bi * bj); + } + } + } + Ok(k) + } + + /// The consistent mass per direction (n × n). + pub fn mass(&self, density: f64, area: f64) -> Result { + let n = self.nodes(); + let mut m = Matrix::zeros(n); + for p in self.points()? { + let factor = density * area * p.weight * p.jacobian; + for i in 0..n { + for j in 0..n { + m.add(i, j, factor * p.n[i] * p.n[j]); + } + } + } + Ok(m) + } + + /// The HRZ lumped mass per node: the consistent diagonal scaled to the member's mass. + pub fn lumped_mass(&self, density: f64, area: f64) -> Result> { + let m = self.mass(density, area)?; + let n = self.nodes(); + let diagonal: Vec = (0..n).map(|i| m.get(i, i)).collect(); + let total = density * area * self.length()?; + let sum: f64 = diagonal.iter().sum(); + if !(sum > 0.0) { + return Err(Error::Degenerate("the member has no mass to lump".into())); + } + Ok(diagonal.iter().map(|d| d * total / sum).collect()) + } + + /// Consistent nodal forces of a force `b` per volume. + pub fn body_force(&self, b: [f64; 3], area: f64) -> Result> { + let mut f = vec![[0.0; 3]; self.nodes()]; + for p in self.points()? { + for (i, fi) in f.iter_mut().enumerate() { + let w = area * p.weight * p.jacobian * p.n[i]; + for k in 0..3 { + fi[k] += w * b[k]; + } + } + } + Ok(f) + } + + /// At each integration point, the row that gives the axial strain from the element's + /// displacements (3n, node by node), and the tangent there: ε = row · u. + pub fn strain_operator(&self) -> Result, [f64; 3])>> { + let n = self.nodes(); + Ok(self.points()?.iter().map(|p| ((0..3 * n).map(|a| p.ds[a / 3] * p.tangent[a % 3]).collect(), p.tangent)).collect()) + } + + /// The axial strain at each integration point for nodal displacements `u`. + pub fn axial_strains(&self, u: &[[f64; 3]]) -> Result> { + if u.len() != self.nodes() { + return Err(Error::Invalid(format!("{} displacements for a {}", u.len(), self.kind.name()))); + } + Ok(self.points()?.iter().map(|p| (0..self.nodes()).map(|i| p.ds[i] * (0..3).map(|k| p.tangent[k] * u[i][k]).sum::()).sum()).collect()) + } + + /// Strain and stress at the integration points as global tensors (see the module docs). + pub fn stresses(&self, young: f64, u: &[[f64; 3]]) -> Result> { + let strains = self.axial_strains(u)?; + Ok(self + .points()? + .iter() + .zip(strains) + .map(|(p, e)| PointStress { xyz: p.xyz, strain: uniaxial(e, p.tangent, 2.0), stress: uniaxial(young * e, p.tangent, 1.0) }) + .collect()) + } + + /// Integration-point values extrapolated to the nodes (the polynomial through the points, + /// evaluated at the nodes' ξ: −1, 1, then 0). + pub fn extrapolate(&self, values: &[[f64; 6]]) -> Result> { + let xs = self.abscissae(); + if values.len() != xs.len() { + return Err(Error::Invalid(format!("{} values for {} integration points", values.len(), xs.len()))); + } + let at = |xi: f64| { + let mut out = [0.0; 6]; + for (a, va) in xs.iter().zip(values) { + let l: f64 = xs.iter().filter(|b| *b != a).map(|b| (xi - b) / (a - b)).product(); + for k in 0..6 { + out[k] += l * va[k]; + } + } + out + }; + Ok([-1.0, 1.0, 0.0].into_iter().take(self.nodes()).map(at).collect()) + } + + /// The tangent at each integration point. + pub fn tangents(&self) -> Result> { + Ok(self.points()?.iter().map(|p| p.tangent).collect()) + } + + /// An upper bound of the highest natural frequency squared of the element with its lumped + /// mass, for explicit time integration: exact (4c²/L²) for rod2, and Gershgorin's bound on + /// M^{-1/2} K M^{-1/2} for rod3. + pub fn frequency_bound(&self, young: f64, density: f64) -> Result { + if self.kind == ElementKind::Rod2 { + let l = chord(self.xyz); + return Ok(4.0 * young / (density * l * l)); + } + // The area cancels: K ∝ A and M ∝ A. + let k = self.stiffness(young, 1.0)?; + let m = self.lumped_mass(density, 1.0)?; + let mut worst: f64 = 0.0; + for i in 0..k.n { + let mi = m.get(i / 3).copied().unwrap_or(0.0); + let row: f64 = (0..k.n).map(|j| k.get(i, j).abs() / (mi * m.get(j / 3).copied().unwrap_or(0.0)).sqrt()).sum(); + worst = worst.max(row); + } + Ok(worst) + } +} + +/// The distance between the ends. +fn chord(xyz: &[[f64; 3]]) -> f64 { + match (xyz.first(), xyz.get(1)) { + (Some(a), Some(b)) => ((a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2) + (a[2] - b[2]).powi(2)).sqrt(), + _ => 0.0, + } +} + +/// The geometry at `xi` of the element on `xyz` (two or three nodes). +fn at(xyz: &[[f64; 3]], xi: f64, weight: f64) -> Result { + let (n, dn) = shape(xyz.len(), xi); + let mut d = [0.0; 3]; + let mut x = [0.0; 3]; + for (i, p) in xyz.iter().enumerate() { + for k in 0..3 { + d[k] += dn[i] * p[k]; + x[k] += n[i] * p[k]; + } + } + let jacobian = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt(); + if !(jacobian > 0.0 && jacobian.is_finite()) { + return Err(Error::Degenerate("the member has no length".into())); + } + Ok(Point { n, ds: dn.map(|v| v / jacobian), tangent: d.map(|c| c / jacobian), jacobian, weight, xyz: x }) +} + +/// The tensor s t tᵀ as a 6-vector (xx, yy, zz, xy, xz, yz), shear terms times `shear` (2 for +/// engineering shear strains). +pub fn uniaxial(s: f64, t: [f64; 3], shear: f64) -> [f64; 6] { + [s * t[0] * t[0], s * t[1] * t[1], s * t[2] * t[2], shear * s * t[0] * t[1], shear * s * t[0] * t[2], shear * s * t[1] * t[2]] +} + +/// The axial value (stress or strain) of a uniaxial tensor along a unit tangent: tᵀ σ t, with +/// engineering shear strains halved when `engineering`. +pub fn axial(v: &[f64; 6], t: [f64; 3], engineering: bool) -> f64 { + let h = if engineering { 0.5 } else { 1.0 }; + v[0] * t[0] * t[0] + v[1] * t[1] * t[1] + v[2] * t[2] * t[2] + 2.0 * h * (v[3] * t[0] * t[1] + v[4] * t[0] * t[2] + v[5] * t[1] * t[2]) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn norm(a: [f64; 3]) -> f64 { + (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt() + } + + /// A rod along a skew direction, and a rod3 with its middle node at the midpoint. + fn straight(kind: ElementKind) -> Vec<[f64; 3]> { + let (a, b) = ([0.5, -1.0, 2.0], [2.5, 1.0, 3.0]); + let mut x = vec![a, b]; + if kind == ElementKind::Rod3 { + x.push([0, 1, 2].map(|k| (a[k] + b[k]) / 2.0)); + } + x + } + + #[test] + fn straight_rods_match_the_textbook_matrices() { + let (young, area, density) = (200.0, 3.0, 2.0); + for kind in [ElementKind::Rod2, ElementKind::Rod3] { + let x = straight(kind); + let t = Truss::new(kind, &x).unwrap(); + let l = norm([2.0, 2.0, 1.0]); + assert!((t.length().unwrap() - l).abs() < 1e-14); + // Rigid translations and the rotation about the axis' normals: zero force. + let k = t.stiffness(young, area).unwrap(); + for d in 0..3 { + let u: Vec = (0..k.n).map(|a| if a % 3 == d { 1.0 } else { 0.0 }).collect(); + assert!(k.mul_vec(&u).iter().all(|f| f.abs() < 1e-12), "{kind:?}"); + } + // A uniform stretch ε along the axis: strain ε at every point, end forces ±EAε. + let eps = 1e-3; + let dir = [2.0 / l, 2.0 / l, 1.0 / l]; + let u: Vec<[f64; 3]> = + x.iter().map(|p| [0, 1, 2].map(|c| eps * dir[c] * ((p[0] - 0.5) * dir[0] + (p[1] + 1.0) * dir[1] + (p[2] - 2.0) * dir[2]))).collect(); + for e in t.axial_strains(&u).unwrap() { + assert!((e - eps).abs() < 1e-15, "{kind:?} {e}"); + } + let flat: Vec = u.iter().flatten().copied().collect(); + let f = k.mul_vec(&flat); + for c in 0..3 { + assert!((f[3 + c] - young * area * eps * dir[c]).abs() < 1e-12, "{kind:?}"); + } + let s = t.stresses(young, &u).unwrap(); + for p in &s { + assert!((axial(&p.stress, dir, false) - young * eps).abs() < 1e-12); + assert!((axial(&p.strain, dir, true) - eps).abs() < 1e-15); + } + let nodal = t.extrapolate(&s.iter().map(|p| p.stress).collect::>()).unwrap(); + assert!(nodal.iter().all(|v| (axial(v, dir, false) - young * eps).abs() < 1e-12)); + // Mass: total ρAL per direction; lumped HRZ. + let m = t.mass(density, area).unwrap(); + let total = density * area * l; + assert!((m.data.iter().sum::() - total).abs() < 1e-12); + let lumped = t.lumped_mass(density, area).unwrap(); + let want: Vec = if kind == ElementKind::Rod2 { vec![0.5, 0.5] } else { vec![1.0 / 6.0, 1.0 / 6.0, 2.0 / 3.0] }; + for (a, b) in lumped.iter().zip(want) { + assert!((a - b * total).abs() < 1e-12, "{kind:?} {lumped:?}"); + } + let fb = t.body_force([0.0, 0.0, -1.0], area).unwrap(); + assert!((fb.iter().map(|f| f[2]).sum::() + area * l).abs() < 1e-12); + // The highest frequency of the free element with its lumped mass: exact for rod2, + // an upper bound of the exact one for rod3. + let bound = t.frequency_bound(young, density).unwrap(); + let c2 = young / density; + if kind == ElementKind::Rod2 { + assert!((bound - 4.0 * c2 / (l * l)).abs() < 1e-9 * bound); + } else { + // Exact: axial eigenproblem of K = EA/(3L)[7 1 −8; 1 7 −8; −8 −8 16] with + // M = ρAL diag(1/6, 1/6, 2/3): ω²max = 24 c²/L² (antisymmetric 12, symmetric 24). + assert!(bound >= 24.0 * c2 / (l * l) * (1.0 - 1e-12), "{bound} vs {}", 24.0 * c2 / (l * l)); + } + } + } + + #[test] + fn bad_geometry_is_an_error() { + assert!(Truss::new(ElementKind::Rod2, &[[1.0; 3], [1.0; 3]]).is_err()); + // The middle node beyond an end folds the member back. + assert!(Truss::new(ElementKind::Rod3, &[[0.0; 3], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]]).is_err()); + assert!(Truss::new(ElementKind::Rod3, &[[0.0; 3], [1.0, 0.0, 0.0]]).is_err()); + assert!(Truss::new(ElementKind::Hex8, &[[0.0; 3]; 8]).is_err()); + assert!(Truss::new(ElementKind::Rod2, &[[0.0; 3], [f64::NAN, 0.0, 0.0]]).is_err()); + assert!(young_modulus(&Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [1.0; 3] }).is_err()); + } + + /// A curved member: a circular arc of radius R (middle node on the arc). Any linear + /// displacement field is reproduced exactly: a uniform expansion u = s x (a ring growing + /// radially by δ = s R) strains it by s everywhere, an infinitesimal rotation not at all. + #[test] + fn curved_rod3_reproduces_linear_fields() { + let (radius, young) = (3.0, 70.0); + for span in [0.3f64, 1.0, 2.0] { + let at = |a: f64| [radius * a.cos(), radius * a.sin(), 0.5]; + let x = [at(0.0), at(span), at(span / 2.0)]; + let t = Truss::new(ElementKind::Rod3, &x).unwrap(); + let (s, w) = (2e-3, [1e-3, -2e-3, 3e-3]); + let u: Vec<[f64; 3]> = + x.iter().map(|p| [0, 1, 2].map(|k| s * p[k] + w[(k + 1) % 3] * p[(k + 2) % 3] - w[(k + 2) % 3] * p[(k + 1) % 3])).collect(); + let points = t.stresses(young, &u).unwrap(); + for (p, tangent) in points.iter().zip(t.tangents().unwrap()) { + assert!((axial(&p.strain, tangent, true) - s).abs() < 1e-15, "{span}"); + assert!((axial(&p.stress, tangent, false) - young * s).abs() < 1e-13, "{span}"); + } + // The arc length: exact for a quadratic arc only as the span shrinks. + let err = (t.length().unwrap() - radius * span).abs() / (radius * span); + assert!(err < 0.01 * span.powi(4), "{span}: {err}"); + } + } +} diff --git a/crates/io/src/mesh.rs b/crates/io/src/mesh.rs index f21a285..8a78849 100644 --- a/crates/io/src/mesh.rs +++ b/crates/io/src/mesh.rs @@ -38,7 +38,7 @@ const PART: &str = "mesh.bin"; const NO_TAG: u32 = u32::MAX; /// The element kind codes of the format. -pub const KIND_CODES: [(ElementKind, u8); 11] = [ +pub const KIND_CODES: [(ElementKind, u8); 12] = [ (ElementKind::Rod2, 1), (ElementKind::Beam2, 2), (ElementKind::Tri3, 3), @@ -50,6 +50,7 @@ pub const KIND_CODES: [(ElementKind, u8); 11] = [ (ElementKind::Wedge6, 9), (ElementKind::Hex8, 10), (ElementKind::Hex20, 11), + (ElementKind::Rod3, 12), ]; fn kind_code(kind: ElementKind) -> u8 { diff --git a/crates/mesh/src/quality/mod.rs b/crates/mesh/src/quality/mod.rs index 18ceecf..efc99e3 100644 --- a/crates/mesh/src/quality/mod.rs +++ b/crates/mesh/src/quality/mod.rs @@ -205,7 +205,9 @@ pub fn element_quality(kind: ElementKind, xyz: &[[f64; 3]]) -> Result return Err(Error::Invalid(format!("quality metrics are not defined for {}", kind.name()))), + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => { + return Err(Error::Invalid(format!("quality metrics are not defined for {}", kind.name()))); + } ElementKind::Tet4 | ElementKind::Tet10 => shape::tet(corners(&p)?), ElementKind::Wedge6 => { let c: [[f64; 3]; 6] = corners(&p)?; diff --git a/crates/model/src/check.rs b/crates/model/src/check.rs index 4146575..76b881b 100644 --- a/crates/model/src/check.rs +++ b/crates/model/src/check.rs @@ -71,6 +71,11 @@ pub struct Resolved { /// apply, so every one of them distributes a load the same way. A force or traction gives /// force components, a moment moment components. pub nodal_loads: BTreeMap>, + /// Every point mass at every node of its target, by definition name, then node. + pub masses: Vec, + /// Every spring and dashpot at every node of its target (or between its two ends), by + /// definition name, then node. + pub springs: Vec, /// The check's warnings (a resolved model has no errors). pub warnings: Vec, /// How each element of a 2D section (plane stress, plane strain, axisymmetric) extends out @@ -97,6 +102,11 @@ impl Resolved { self.sets.get(set).map(|s| s.edges.iter().filter(|e| self.planar.contains_key(&e.element)).copied().collect()).unwrap_or_default() } + /// The node labels point masses, springs and dashpots act on (repeats possible). + pub fn discrete_labels(&self) -> impl Iterator + '_ { + self.masses.iter().map(|m| m.node).chain(self.springs.iter().flat_map(|s| std::iter::once(s.first).chain(s.second))) + } + /// The nodes a boundary condition, force or moment acts on. pub fn target_nodes(&self, target: &Target) -> &[u32] { match target { @@ -160,9 +170,10 @@ impl Model { let mut resolved = Resolved::default(); let units = self.units; - // A model of fluid steps alone needs no mesh: its domains carry their own grid. + // A model of fluid steps alone needs no mesh: its domains carry their own grid; one of + // point masses, springs and dashpots alone is its own mesh. let fluid_only = !self.steps.is_empty() && self.steps.iter().all(|s| s.procedure.fluid_domain().is_some()); - if self.elements.is_empty() && !fluid_only { + if self.elements.is_empty() && !fluid_only && self.masses.is_empty() && self.springs.is_empty() { d.error("model", "the mesh has no elements"); } if units.length().is_none() { @@ -229,7 +240,7 @@ impl Model { k.dimension() == 2 } SectionKind::Beam { .. } => k == ElementKind::Beam2, - SectionKind::Rod { .. } => k == ElementKind::Rod2, + SectionKind::Rod { .. } => matches!(k, ElementKind::Rod2 | ElementKind::Rod3), }; // The linear element with an incompatible-modes formulation for this kind of section. let incompatible_kind = if sec.kind.is_planar() { ElementKind::Quad4 } else { ElementKind::Hex8 }; @@ -362,6 +373,9 @@ impl Model { } } + // Point masses, springs and dashpots. + let (masses, springs) = self.resolve_discrete(&resolved, &mut d, &has_rotations); + // Loads. let mut nodal_loads = BTreeMap::new(); for (name, load) in &self.loads { @@ -449,11 +463,123 @@ impl Model { } resolved.nodal_loads = nodal_loads; + resolved.masses = masses; + resolved.springs = springs; let mut out = d.0; out.sort_by_key(|a| a.severity); (resolved, out) } + /// Point masses, springs and dashpots expanded to their nodes, reporting what cannot be + /// expanded. + fn resolve_discrete( + &self, + resolved: &Resolved, + d: &mut Diagnostics, + has_rotations: &dyn Fn(&Target) -> bool, + ) -> (Vec, Vec) { + use crate::{Direction, ResolvedDirection, ResolvedMass, ResolvedSpring}; + let mut masses = Vec::new(); + let mut springs = Vec::new(); + let mut count = 0usize; + let mut over = false; + let mut take = |n: usize| { + count = count.saturating_add(n); + over |= count > crate::discrete::MAX_DISCRETE; + !over + }; + for (name, m) in &self.masses { + let subject = format!("point mass `{name}`"); + let nodes = resolved.target_nodes(&m.target); + if nodes.is_empty() { + d.error(&subject, format!("{} has no nodes", m.target)); + continue; + } + if m.inertia.is_some_and(|i| i.iter().any(|c| *c != 0.0)) && !has_rotations(&m.target) { + d.error(&subject, format!("{} has no nodes with rotations for its inertia; put it on a reference point", m.target)); + } + if !take(nodes.len()) { + break; + } + let mass = m.mass.map_or([0.0; 3], crate::Mass::components); + let inertia = m.inertia.unwrap_or_default(); + masses.extend(nodes.iter().map(|n| ResolvedMass { name: name.clone(), node: *n, mass, inertia })); + } + for (name, s) in &self.springs { + let subject = format!("{} `{name}`", s.kind()); + let first = resolved.target_nodes(s.first()); + let rotational = s.direction().is_rotational(); + if first.is_empty() { + d.error(&subject, format!("{} has no nodes", s.first())); + continue; + } + if rotational && s.targets().any(|t| !has_rotations(t)) { + d.error(&subject, "acts on a rotation, and its nodes have none; connect reference points (or shell and beam nodes)"); + } + let unit = |v: [f64; 3]| { + let len = (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt(); + (len > 0.0 && len.is_finite()).then(|| v.map(|c| c / len)) + }; + let pairs: Vec<(u32, Option)> = match s.second() { + None => first.iter().map(|n| (*n, None)).collect(), + Some(t) => { + let second = resolved.target_nodes(t); + if first.len() != 1 || second.len() != 1 { + d.error( + &subject, + format!( + "a two-point {} connects one node to one node; {} has {} and {t} has {} (use a vertex set or a reference point at each end)", + s.kind(), + s.first(), + first.len(), + second.len() + ), + ); + continue; + } + match (first.first(), second.first()) { + (Some(a), Some(b)) if a == b => { + d.error(&subject, format!("both ends are node {a}")); + continue; + } + (Some(a), Some(b)) => vec![(*a, Some(*b))], + _ => continue, + } + } + }; + if !take(pairs.len()) { + break; + } + for (a, b) in pairs { + let direction = match s.direction() { + Direction::Dof(k) => ResolvedDirection::Dof(k), + Direction::Vector(v) => match unit(v) { + Some(n) => ResolvedDirection::Vector(n), + None => { + d.error(&subject, "its direction vector has no length"); + break; + } + }, + Direction::Line => { + let ends = b.and_then(|b| Some((self.label_xyz(a)?, self.label_xyz(b)?))); + match ends.and_then(|(x, y)| unit([y[0] - x[0], y[1] - x[1], y[2] - x[2]])) { + Some(n) => ResolvedDirection::Line(n), + None => { + d.error(&subject, "its two ends are at the same point, so it has no line to act along"); + break; + } + } + } + }; + springs.push(ResolvedSpring { name: name.clone(), dashpot: s.is_dashpot(), first: a, second: b, direction, value: s.value() }); + } + } + if over { + d.error("model", format!("point masses, springs and dashpots resolve to more than {} copies", crate::discrete::MAX_DISCRETE)); + } + (masses, springs) + } + /// A force, moment or traction as nodal values (`None` for other loads), or why it cannot be /// distributed. fn nodal_values(&self, load: &Load, r: &Resolved) -> std::result::Result>, String> { @@ -505,13 +631,17 @@ impl Model { })) } - /// The rules of a 2D model: one kind of 2D model (plane or axisymmetric) and nothing else, in - /// the XY plane (x ≥ 0 when axisymmetric), elements numbered counter-clockwise. + /// The rules of a 2D model: one kind of 2D model (plane or axisymmetric) and nothing else but + /// trusses (plane models only), in the XY plane (x ≥ 0 when axisymmetric), elements numbered + /// counter-clockwise. fn check_planar(&self, r: &Resolved, d: &mut Diagnostics) { if !r.is_planar() { return; } - let others: Vec<&str> = self.sections.iter().filter(|(_, s)| !s.kind.is_planar()).map(|(n, _)| n.as_str()).collect(); + let axisymmetric = r.extent() == Some(Extent::Revolved); + // Trusses join a plane model with its two in-plane unknowns per node. + let truss = |s: &crate::Section| matches!(s.kind, SectionKind::Rod { .. }) && !axisymmetric; + let others: Vec<&str> = self.sections.iter().filter(|(_, s)| !s.kind.is_planar() && !truss(s)).map(|(n, _)| n.as_str()).collect(); if !others.is_empty() { d.error( "model", @@ -521,7 +651,6 @@ impl Model { ), ); } - let axisymmetric = r.extent() == Some(Extent::Revolved); if axisymmetric && r.planar.values().any(|e| *e != Extent::Revolved) { d.error("model", "axisymmetric and plane (stress or strain) sections cannot be mixed in one model"); } @@ -542,6 +671,14 @@ impl Model { clockwise.push(label); } } + for sec in self.sections.values().filter(|s| truss(s)) { + for label in r.sets.get(&sec.set).map(|s| s.elements.as_slice()).unwrap_or_default() { + let Some(e) = self.elements.get(label) else { continue }; + if e.nodes.iter().filter_map(|n| self.nodes.get(n)).any(|n| n.xyz[2] != 0.0) { + off_plane.push(*label); + } + } + } if !off_plane.is_empty() { d.error("model", format!("2D elements {} have nodes off the XY plane (z ≠ 0); 2D models lie in the XY plane", first_labels(&off_plane))); } diff --git a/crates/model/src/constraints.rs b/crates/model/src/constraints.rs index 764d2b2..e0bda58 100644 --- a/crates/model/src/constraints.rs +++ b/crates/model/src/constraints.rs @@ -104,6 +104,7 @@ impl Model { .map(|(n, _)| format!("boundary condition `{n}`")), ); users.extend(self.loads.iter().filter(|(_, l)| l.reference_point() == Some(name)).map(|(n, _)| format!("load `{n}`"))); + users.extend(self.discrete_users(|t| t.reference_point() == Some(name))); not_in_use("reference point", name, users)?; remove_existing(&mut self.reference_points, "reference point", name).map(|_| ()) } diff --git a/crates/model/src/discrete.rs b/crates/model/src/discrete.rs new file mode 100644 index 0000000..a254b19 --- /dev/null +++ b/crates/model/src/discrete.rs @@ -0,0 +1,501 @@ +//! Discrete elements: point masses with rotary inertia, springs and dashpots. +//! +//! They are definitions, not mesh elements: each acts on a target (the nodes of a set, or a +//! reference point) and survives remeshing like a load does. Every node of the target gets its +//! own copy, with the full value: a point mass of 2 on a set of five nodes adds 2 at each of them +//! (a mass shared over a face is a point mass on a reference point coupled to the face). +//! +//! - A **point mass** adds a translational mass, the same in every direction or one per global +//! direction, and a rotary inertia tensor about the node (which needs rotational degrees of +//! freedom: a reference point, or a shell or beam node). +//! - A **spring** resists the relative displacement of its two ends along a direction with a +//! force `stiffness × relative displacement`; a **dashpot** resists the relative velocity with +//! `coefficient × relative velocity`. With no second end, the spring or dashpot connects each +//! node of the target to the ground. The direction is a global degree of freedom (1–3 +//! translations, 4–6 rotations: then a torsional spring or dashpot), a fixed vector, or the +//! line through the two ends (two-point only, from the undeformed positions). +//! +//! The linear (small displacement) form is the textbook spring element: for a unit direction n, +//! the element matrix is k [n nᵀ, −n nᵀ; −n nᵀ, n nᵀ] on the two ends' translations (R. D. Cook, +//! D. S. Malkus, M. E. Plesha, R. J. Witt, *Concepts and Applications of Finite Element +//! Analysis*, 4th ed., Wiley 2002, §2.2 and §2.10), and k n nᵀ for a grounded one. Dashpots +//! have the same form on velocities; static and frequency steps ignore them. + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::bcs::Target; +use crate::{Error, Model, Result, context, finite, insert_new, remove_existing, replace_existing}; + +/// Most copies of discrete elements a model resolves to (definitions × target nodes). +pub const MAX_DISCRETE: usize = 10_000_000; + +/// A translational mass: the same in every direction, or one per global direction. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(untagged)] +pub enum Mass { + Isotropic(f64), + Directional([f64; 3]), +} + +impl Mass { + /// The mass in X, Y and Z. + pub fn components(self) -> [f64; 3] { + match self { + Mass::Isotropic(m) => [m; 3], + Mass::Directional(m) => m, + } + } +} + +/// A point mass and rotary inertia at every node of a target. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct PointMass { + pub target: Target, + #[serde(default, skip_serializing_if = "Option::is_none")] + pub mass: Option, + /// The inertia tensor about the node in global axes, `[I11, I22, I33, I12, I13, I23]`: the + /// components of the symmetric tensor J with angular momentum H = J ω (products of inertia + /// carry their sign: I12 = −∫ x y dm). + #[serde(default, skip_serializing_if = "Option::is_none")] + pub inertia: Option<[f64; 6]>, +} + +impl PointMass { + pub fn validate(&self) -> Result<()> { + if self.mass.is_none() && self.inertia.is_none() { + return Err(Error::Invalid("a point mass needs a mass, an inertia or both".into())); + } + if let Some(m) = self.mass { + for c in m.components() { + if finite("mass", c)? < 0.0 { + return Err(Error::Invalid(format!("a mass is not negative, got {c}"))); + } + } + if m.components().iter().all(|c| *c == 0.0) && self.inertia.is_none() { + return Err(Error::Invalid("a point mass of zero adds nothing".into())); + } + } + if let Some(i) = self.inertia { + for c in i { + finite("inertia", c)?; + } + // Positive semidefinite: every principal minor is non-negative (to round-off). + let j = inertia_matrix(i); + let scale = i.iter().fold(0.0f64, |m, c| m.max(c.abs())); + let tolerance = 1e-12 * scale; + let minors = [ + j[0][0], + j[1][1], + j[2][2], + j[0][0] * j[1][1] - j[0][1] * j[0][1], + j[0][0] * j[2][2] - j[0][2] * j[0][2], + j[1][1] * j[2][2] - j[1][2] * j[1][2], + j[0][0] * (j[1][1] * j[2][2] - j[1][2] * j[2][1]) - j[0][1] * (j[1][0] * j[2][2] - j[1][2] * j[2][0]) + + j[0][2] * (j[1][0] * j[2][1] - j[1][1] * j[2][0]), + ]; + let powers = [1, 1, 1, 2, 2, 2, 3]; + if minors.iter().zip(powers).any(|(m, p)| *m < -tolerance * scale.powi(p - 1)) { + return Err(Error::Invalid( + "an inertia tensor must be positive semidefinite (non-negative moments, products no larger than the moments allow)".into(), + )); + } + if i.iter().all(|c| *c == 0.0) && self.mass.is_none_or(|m| m.components().iter().all(|c| *c == 0.0)) { + return Err(Error::Invalid("a point mass of zero adds nothing".into())); + } + } + Ok(()) + } +} + +/// `[I11, I22, I33, I12, I13, I23]` as a symmetric matrix. +pub fn inertia_matrix(i: [f64; 6]) -> [[f64; 3]; 3] { + let [i11, i22, i33, i12, i13, i23] = i; + [[i11, i12, i13], [i12, i22, i23], [i13, i23, i33]] +} + +/// The direction a spring or dashpot acts in. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum Direction { + /// A global degree of freedom: 1–3 translations along X, Y, Z, 4–6 rotations about them. + Dof(u8), + /// A fixed direction in global axes (translations; normalized when used). + Vector([f64; 3]), + /// The line through the two ends' undeformed positions (two-point only). + Line, +} + +impl Direction { + /// Whether it acts on rotations. + pub fn is_rotational(self) -> bool { + matches!(self, Direction::Dof(4..=6)) + } + + fn validate(self, two_point: bool) -> Result<()> { + match self { + Direction::Dof(d) if !(1..=6).contains(&d) => Err(Error::Invalid(format!("degree of freedom {d} is not in 1..=6"))), + Direction::Dof(_) => Ok(()), + Direction::Vector(v) => { + for c in v { + finite("direction", c)?; + } + if v.iter().all(|c| *c == 0.0) { + return Err(Error::Invalid("a direction vector must not be zero".into())); + } + Ok(()) + } + Direction::Line if !two_point => Err(Error::Invalid("the direction `line` needs a second end (it runs between the two ends)".into())), + Direction::Line => Ok(()), + } + } +} + +/// A linear spring or dashpot between the nodes of `first` and the ground, or between `first` +/// and `second` (each one node: a vertex set or a reference point). +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] +pub enum Spring { + /// Force (or moment) = stiffness × relative displacement (or rotation). + Spring { + first: Target, + #[serde(default, skip_serializing_if = "Option::is_none")] + second: Option, + direction: Direction, + stiffness: f64, + }, + /// Force (or moment) = coefficient × relative velocity (or angular velocity); dynamic steps + /// only. + Dashpot { + first: Target, + #[serde(default, skip_serializing_if = "Option::is_none")] + second: Option, + direction: Direction, + coefficient: f64, + }, +} + +impl Spring { + pub fn first(&self) -> &Target { + match self { + Spring::Spring { first, .. } | Spring::Dashpot { first, .. } => first, + } + } + + pub fn second(&self) -> Option<&Target> { + match self { + Spring::Spring { second, .. } | Spring::Dashpot { second, .. } => second.as_ref(), + } + } + + pub fn direction(&self) -> Direction { + match self { + Spring::Spring { direction, .. } | Spring::Dashpot { direction, .. } => *direction, + } + } + + /// The stiffness or the damping coefficient. + pub fn value(&self) -> f64 { + match self { + Spring::Spring { stiffness, .. } => *stiffness, + Spring::Dashpot { coefficient, .. } => *coefficient, + } + } + + pub fn is_dashpot(&self) -> bool { + matches!(self, Spring::Dashpot { .. }) + } + + /// `spring` or `dashpot`. + pub fn kind(&self) -> &'static str { + if self.is_dashpot() { "dashpot" } else { "spring" } + } + + /// The targets it acts on. + pub fn targets(&self) -> impl Iterator { + std::iter::once(self.first()).chain(self.second()) + } + + pub fn validate(&self) -> Result<()> { + let what = if self.is_dashpot() { "damping coefficient" } else { "stiffness" }; + let v = finite(what, self.value())?; + if v <= 0.0 { + return Err(Error::Invalid(format!("a {} needs a positive {what}, got {v}", self.kind()))); + } + self.direction().validate(self.second().is_some())?; + if self.second() == Some(self.first()) { + return Err(Error::Invalid(format!("a {} cannot connect a target to itself", self.kind()))); + } + Ok(()) + } +} + +/// A resolved direction: a global degree of freedom (1–6), or a unit vector acting on the +/// translations (from a fixed vector, or along the line of the ends). +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum ResolvedDirection { + Dof(u8), + Vector([f64; 3]), + Line([f64; 3]), +} + +impl ResolvedDirection { + /// The coefficients of the relative motion on each end's six degrees of freedom: the + /// spring's elongation is Σ c (u_second − u_first). + pub fn coefficients(self) -> [f64; 6] { + let mut c = [0.0; 6]; + match self { + ResolvedDirection::Dof(d) => { + if let Some(x) = usize::from(d).checked_sub(1).and_then(|i| c.get_mut(i)) { + *x = 1.0; + } + } + ResolvedDirection::Vector(n) | ResolvedDirection::Line(n) => c[..3].copy_from_slice(&n), + } + c + } + + pub fn is_rotational(self) -> bool { + matches!(self, ResolvedDirection::Dof(4..=6)) + } + + /// The unit vector it acts along (for a rotation, the axis it turns about). + pub fn axis(self) -> [f64; 3] { + let c = self.coefficients(); + if self.is_rotational() { [c[3], c[4], c[5]] } else { [c[0], c[1], c[2]] } + } +} + +/// One point mass at one node. +#[derive(Debug, Clone, PartialEq)] +pub struct ResolvedMass { + /// The definition it comes from. + pub name: String, + pub node: u32, + /// Mass in X, Y, Z. + pub mass: [f64; 3], + /// `[I11, I22, I33, I12, I13, I23]` (zero when there is none). + pub inertia: [f64; 6], +} + +impl ResolvedSpring { + /// `spring` or `dashpot`. + pub fn kind(&self) -> &'static str { + if self.dashpot { "dashpot" } else { "spring" } + } +} + +/// One spring or dashpot at one node (grounded) or between two nodes. +#[derive(Debug, Clone, PartialEq)] +pub struct ResolvedSpring { + pub name: String, + pub dashpot: bool, + pub first: u32, + /// `None`: to the ground. + pub second: Option, + pub direction: ResolvedDirection, + /// Stiffness, or damping coefficient. + pub value: f64, +} + +impl Model { + pub fn masses(&self) -> &BTreeMap { + &self.masses + } + + pub fn mass(&self, name: &str) -> Option<&PointMass> { + self.masses.get(name) + } + + pub fn springs(&self) -> &BTreeMap { + &self.springs + } + + pub fn spring(&self, name: &str) -> Option<&Spring> { + self.springs.get(name) + } + + fn check_mass(&self, name: &str, m: &PointMass) -> Result<()> { + m.validate().map_err(|e| context("point mass", name, e))?; + m.target.check_exists(self).map_err(|e| context("point mass", name, e)) + } + + pub fn add_mass(&mut self, name: &str, mass: PointMass) -> Result<()> { + self.check_mass(name, &mass)?; + insert_new(&mut self.masses, "point mass", name, mass) + } + + pub fn replace_mass(&mut self, name: &str, mass: PointMass) -> Result<()> { + self.check_mass(name, &mass)?; + replace_existing(&mut self.masses, "point mass", name, mass) + } + + pub fn remove_mass(&mut self, name: &str) -> Result<()> { + remove_existing(&mut self.masses, "point mass", name).map(|_| ()) + } + + fn check_spring(&self, name: &str, s: &Spring) -> Result<()> { + let what = s.kind(); + s.validate().map_err(|e| context(what, name, e))?; + for t in s.targets() { + t.check_exists(self).map_err(|e| context(what, name, e))?; + } + Ok(()) + } + + pub fn add_spring(&mut self, name: &str, spring: Spring) -> Result<()> { + self.check_spring(name, &spring)?; + insert_new(&mut self.springs, "spring", name, spring) + } + + pub fn replace_spring(&mut self, name: &str, spring: Spring) -> Result<()> { + self.check_spring(name, &spring)?; + replace_existing(&mut self.springs, "spring", name, spring) + } + + pub fn remove_spring(&mut self, name: &str) -> Result<()> { + remove_existing(&mut self.springs, "spring", name).map(|_| ()) + } + + /// The definitions among point masses, springs and dashpots that act on a target. + pub(crate) fn discrete_users(&self, uses: impl Fn(&Target) -> bool) -> Vec { + let mut users: Vec = self.masses.iter().filter(|(_, m)| uses(&m.target)).map(|(n, _)| format!("point mass `{n}`")).collect(); + users.extend(self.springs.iter().filter(|(_, s)| s.targets().any(&uses)).map(|(n, s)| format!("{} `{n}`", s.kind()))); + users + } + + /// The coordinates of a mesh node or a reference point's node label. + pub fn label_xyz(&self, label: u32) -> Option<[f64; 3]> { + self.nodes.get(&label).map(|n| n.xyz).or_else(|| self.reference_points.values().find(|p| p.label == label).map(|p| p.xyz)) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::Set; + + #[test] + fn definitions_validate_and_round_trip() { + let mut m = Model::default(); + m.add_set("tip", Set { nodes: vec![1], ..Set::default() }).unwrap(); + let json = serde_json::json!({"target": {"set": "tip"}, "mass": [1.0, 2.0, 0.0], "inertia": [1.0, 2.0, 3.0, 0.5, 0.0, 0.0]}); + let pm: PointMass = serde_json::from_value(json.clone()).unwrap(); + assert_eq!(pm.mass, Some(Mass::Directional([1.0, 2.0, 0.0]))); + assert_eq!(serde_json::to_value(&pm).unwrap(), json); + m.add_mass("pm", pm).unwrap(); + let iso: PointMass = serde_json::from_value(serde_json::json!({"target": {"set": "tip"}, "mass": 3})).unwrap(); + assert_eq!(iso.mass.map(Mass::components), Some([3.0; 3])); + for bad in [ + serde_json::json!({"target": {"set": "tip"}}), + serde_json::json!({"target": {"set": "tip"}, "mass": -1}), + serde_json::json!({"target": {"set": "tip"}, "mass": 0}), + serde_json::json!({"target": {"set": "tip"}, "inertia": [1, 1, 1, 2, 0, 0]}), + serde_json::json!({"target": {"set": "tip"}, "inertia": [-1, 1, 1, 0, 0, 0]}), + serde_json::json!({"target": {"set": "tip"}, "inertia": [0, 0, 0, 0, 0, 0]}), + ] { + let pm: PointMass = serde_json::from_value(bad.clone()).unwrap(); + assert!(m.add_mass("x", pm).is_err(), "{bad}"); + } + assert!(serde_json::from_value::(serde_json::json!({"target": {"set": "tip"}, "mass": 1, "masss": 2})).is_err()); + let missing: PointMass = serde_json::from_value(serde_json::json!({"target": {"set": "nope"}, "mass": 1})).unwrap(); + assert!(matches!(m.add_mass("x", missing), Err(Error::NotFound(_)))); + + let s: Spring = + serde_json::from_value(serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": {"dof": 3}, "stiffness": 100})) + .unwrap(); + m.add_spring("k", s).unwrap(); + let d: Spring = serde_json::from_value( + serde_json::json!({"type": "dashpot", "first": {"set": "tip"}, "direction": {"vector": [1, 1, 0]}, "coefficient": 2}), + ) + .unwrap(); + assert!(d.is_dashpot()); + m.add_spring("c", d).unwrap(); + for bad in [ + serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": "line", "stiffness": 1}), + serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": {"dof": 7}, "stiffness": 1}), + serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": {"vector": [0, 0, 0]}, "stiffness": 1}), + serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": {"dof": 1}, "stiffness": 0}), + serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": {"dof": 1}, "stiffness": f64::NAN}), + serde_json::json!({"type": "dashpot", "first": {"set": "tip"}, "second": {"set": "tip"}, "direction": {"dof": 1}, "coefficient": 1}), + ] { + let Ok(s) = serde_json::from_value::(bad.clone()) else { continue }; + assert!(m.add_spring("x", s).is_err(), "{bad}"); + } + assert!( + serde_json::from_value::(serde_json::json!({"type": "dashpot", "first": {"set": "tip"}, "direction": "line", "stiffness": 1})) + .is_err() + ); + // The set is in use. + assert!(matches!(m.remove_set("tip"), Err(Error::InUse(_)))); + m.remove_mass("pm").unwrap(); + m.remove_spring("k").unwrap(); + m.remove_spring("c").unwrap(); + m.remove_set("tip").unwrap(); + } + + #[test] + fn directions_give_unit_coefficients() { + assert_eq!(ResolvedDirection::Dof(5).coefficients(), [0.0, 0.0, 0.0, 0.0, 1.0, 0.0]); + assert_eq!(ResolvedDirection::Line([0.6, 0.8, 0.0]).coefficients(), [0.6, 0.8, 0.0, 0.0, 0.0, 0.0]); + assert!(ResolvedDirection::Dof(4).is_rotational() && !ResolvedDirection::Dof(3).is_rotational()); + } + + #[test] + fn resolution_expands_and_checks_targets() { + use crate::{Procedure, ReferencePoint, Step}; + let mut m = Model::default(); + m.set_node(1, [0.0; 3]).unwrap(); + m.set_node(2, [3.0, 4.0, 0.0]).unwrap(); + m.add_set("a", Set { nodes: vec![1], ..Set::default() }).unwrap(); + m.add_set("b", Set { nodes: vec![2], ..Set::default() }).unwrap(); + m.add_set("ab", Set { nodes: vec![1, 2], ..Set::default() }).unwrap(); + m.add_reference_point("RP", ReferencePoint { xyz: [0.0, 0.0, 1.0], label: 0 }).unwrap(); + m.add_step(Step::new("Modes", Procedure::Frequency { modes: 1, min_hz: None, max_hz: None }), None).unwrap(); + let spring = |first: &str, second: Option<&str>, direction| Spring::Spring { + first: Target::Set(first.into()), + second: second.map(|s| Target::Set(s.into())), + direction, + stiffness: 10.0, + }; + m.add_spring("line", spring("a", Some("b"), Direction::Line)).unwrap(); + m.add_spring("ground", spring("ab", None, Direction::Vector([0.0, 0.0, 2.0]))).unwrap(); + m.add_mass("m", PointMass { target: Target::Set("ab".into()), mass: Some(Mass::Isotropic(2.0)), inertia: None }).unwrap(); + m.add_mass("j", PointMass { target: Target::ReferencePoint("RP".into()), mass: None, inertia: Some([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) }) + .unwrap(); + // No mesh elements: the discrete elements are the model. + let r = m.resolve().unwrap(); + assert_eq!(r.springs.len(), 3); + assert_eq!(r.springs[0].direction, ResolvedDirection::Vector([0.0, 0.0, 1.0])); + assert_eq!((r.springs[2].first, r.springs[2].second), (1, Some(2))); + assert_eq!(r.springs[2].direction, ResolvedDirection::Line([0.6, 0.8, 0.0])); + assert_eq!(r.masses.len(), 3); + assert_eq!(r.masses.iter().map(|x| x.node).collect::>(), vec![99_999_999, 1, 2]); + + // Errors: two-point ends with several nodes, rotations at solid nodes, inertia at mesh nodes. + m.add_spring("wide", spring("ab", Some("b"), Direction::Dof(1))).unwrap(); + m.add_spring("twist", spring("a", None, Direction::Dof(6))).unwrap(); + m.add_mass("bad", PointMass { target: Target::Set("a".into()), mass: None, inertia: Some([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) }).unwrap(); + let errors = m.resolve().unwrap_err().to_string(); + for want in ["spring `wide`: a two-point spring connects one node to one node", "spring `twist`: acts on a rotation", "point mass `bad`"] { + assert!(errors.contains(want), "{errors}"); + } + // A reference point used by a discrete element cannot be deleted. + assert!(matches!(m.remove_reference_point("RP"), Err(Error::InUse(_)))); + // Torsion at a reference point is fine. + m.remove_spring("twist").unwrap(); + m.remove_spring("wide").unwrap(); + m.remove_mass("bad").unwrap(); + m.add_spring( + "twist", + Spring::Spring { first: Target::ReferencePoint("RP".into()), second: None, direction: Direction::Dof(6), stiffness: 1.0 }, + ) + .unwrap(); + m.resolve().unwrap(); + // The definitions survive a file round trip. + let back: Model = serde_json::from_str(&serde_json::to_string(&m).unwrap()).unwrap(); + assert_eq!(back, m); + } +} diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index fa64148..5dec6cc 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -8,7 +8,8 @@ //! - named sets ([`sets`]) defined by geometry names or explicit labels, resolved to nodes, //! elements and element faces against the current mesh, so they survive remeshing; //! - materials and the unit system ([`materials`], [`units`]), sections with their element -//! formulation ([`sections`]), reference points and couplings ([`constraints`]); +//! formulation ([`sections`]), reference points and couplings ([`constraints`]), point masses, +//! springs and dashpots ([`discrete`]); //! - boundary conditions ([`bcs`]), loads ([`loads`]) and steps with an explicit table of which //! boundary conditions and loads are active in each step ([`steps`]); //! - the check that names every missing or bad item, and the resolved view the deck writers and @@ -33,6 +34,7 @@ pub mod bcs; pub mod check; pub mod constraints; +pub mod discrete; pub mod fluids; pub mod geometry; pub mod loads; @@ -51,6 +53,7 @@ use serde::{Deserialize, Serialize}; pub use bcs::{BoundaryCondition, Dofs, Symmetry, Target}; pub use check::{Diagnostic, Resolved, ResolvedSet, Severity}; pub use constraints::{Coupling, CouplingKind, ReferencePoint}; +pub use discrete::{Direction, Mass, PointMass, ResolvedDirection, ResolvedMass, ResolvedSpring, Spring}; pub use fluids::{BoxFace, BoxFaces, FaceCondition, FluidBoundary, FluidDomain, Profile, SurfaceBoundary, SurfaceCondition}; pub use geometry::{Classification, Entity, GeometryRef}; pub use loads::{Distribution, Load, centrifugal_acceleration}; @@ -113,6 +116,10 @@ pub struct Model { sections: BTreeMap, reference_points: BTreeMap, couplings: BTreeMap, + #[serde(skip_serializing_if = "BTreeMap::is_empty")] + masses: BTreeMap, + #[serde(skip_serializing_if = "BTreeMap::is_empty")] + springs: BTreeMap, boundary_conditions: BTreeMap, loads: BTreeMap, #[serde(skip_serializing_if = "BTreeMap::is_empty")] @@ -144,6 +151,10 @@ struct ModelData { #[serde(default)] couplings: BTreeMap, #[serde(default)] + masses: BTreeMap, + #[serde(default)] + springs: BTreeMap, + #[serde(default)] boundary_conditions: BTreeMap, #[serde(default)] loads: BTreeMap, @@ -198,6 +209,12 @@ impl TryFrom for Model { for (n, v) in d.couplings { m.add_coupling(&n, v)?; } + for (n, v) in d.masses { + m.add_mass(&n, v)?; + } + for (n, v) in d.springs { + m.add_spring(&n, v)?; + } for (n, v) in d.boundary_conditions { m.add_boundary_condition(&n, v)?; } @@ -307,6 +324,8 @@ impl Model { sections: self.sections.clone(), reference_points: self.reference_points.clone(), couplings: self.couplings.clone(), + masses: self.masses.clone(), + springs: self.springs.clone(), boundary_conditions: self.boundary_conditions.clone(), loads: self.loads.clone(), fluid_domains: self.fluid_domains.clone(), @@ -315,13 +334,15 @@ impl Model { } /// Counts of each kind of definition, for summaries. - pub fn definition_counts(&self) -> [(&'static str, usize); 9] { + pub fn definition_counts(&self) -> [(&'static str, usize); 11] { [ ("sets", self.sets.len()), ("materials", self.materials.len()), ("sections", self.sections.len()), ("reference_points", self.reference_points.len()), ("couplings", self.couplings.len()), + ("masses", self.masses.len()), + ("springs", self.springs.len()), ("boundary_conditions", self.boundary_conditions.len()), ("loads", self.loads.len()), ("fluid_domains", self.fluid_domains.len()), diff --git a/crates/model/src/mesh.rs b/crates/model/src/mesh.rs index 019ae13..c1552cc 100644 --- a/crates/model/src/mesh.rs +++ b/crates/model/src/mesh.rs @@ -14,6 +14,7 @@ pub struct Node { /// Element topology. Node order is corners first, then mid-side nodes, with positive volume for /// a right-handed ordering: +/// - rod3: the two ends, then the middle node; /// - tet10: mid-side nodes 5–10 on edges 1-2, 2-3, 3-1, 1-4, 2-4, 3-4; /// - wedge6: triangle 1-2-3, then 4-5-6 on the side its right-hand normal points to; /// - hex8/hex20: face 1-2-3-4, then 5-6-7-8 above it; hex20's mid-side nodes 9–12 lie on the @@ -25,8 +26,10 @@ pub struct Node { #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] #[serde(rename_all = "snake_case")] pub enum ElementKind { - /// Two-node axial member. + /// Two-node axial member (truss). Rod2, + /// Three-node axial member (quadratic truss, straight or curved). + Rod3, /// Two-node beam. Beam2, /// Three-node triangle (shell or membrane). @@ -50,8 +53,9 @@ pub enum ElementKind { } impl ElementKind { - pub const ALL: [ElementKind; 11] = [ + pub const ALL: [ElementKind; 12] = [ ElementKind::Rod2, + ElementKind::Rod3, ElementKind::Beam2, ElementKind::Tri3, ElementKind::Tri6, @@ -67,7 +71,7 @@ impl ElementKind { pub fn node_count(self) -> usize { match self { ElementKind::Rod2 | ElementKind::Beam2 => 2, - ElementKind::Tri3 => 3, + ElementKind::Rod3 | ElementKind::Tri3 => 3, ElementKind::Quad4 | ElementKind::Tet4 => 4, ElementKind::Tri6 | ElementKind::Wedge6 => 6, ElementKind::Quad8 | ElementKind::Hex8 => 8, @@ -79,7 +83,7 @@ impl ElementKind { /// Corner nodes: the first ones in the node list (the rest are mid-side nodes). pub fn corner_count(self) -> usize { match self { - ElementKind::Rod2 | ElementKind::Beam2 => 2, + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => 2, ElementKind::Tri3 | ElementKind::Tri6 => 3, ElementKind::Quad4 | ElementKind::Quad8 | ElementKind::Tet4 | ElementKind::Tet10 => 4, ElementKind::Wedge6 => 6, @@ -95,7 +99,7 @@ impl ElementKind { /// Topological dimension: 1 for line elements, 2 for shells, 3 for solids. pub fn dimension(self) -> u8 { match self { - ElementKind::Rod2 | ElementKind::Beam2 => 1, + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => 1, ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => 2, ElementKind::Tet4 | ElementKind::Tet10 | ElementKind::Wedge6 | ElementKind::Hex8 | ElementKind::Hex20 => 3, } @@ -105,6 +109,7 @@ impl ElementKind { pub fn name(self) -> &'static str { match self { ElementKind::Rod2 => "rod2", + ElementKind::Rod3 => "rod3", ElementKind::Beam2 => "beam2", ElementKind::Tri3 => "tri3", ElementKind::Tri6 => "tri6", @@ -134,7 +139,7 @@ impl ElementKind { /// - Line elements have none. pub fn faces(self) -> &'static [&'static [usize]] { match self { - ElementKind::Rod2 | ElementKind::Beam2 => &[], + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => &[], ElementKind::Tri3 => &[&[0, 1, 2], &[0, 2, 1]], ElementKind::Tri6 => &[&[0, 1, 2, 3, 4, 5], &[0, 2, 1, 5, 4, 3]], ElementKind::Quad4 => &[&[0, 1, 2, 3], &[0, 3, 2, 1]], @@ -160,6 +165,7 @@ impl ElementKind { pub fn edges(self) -> &'static [&'static [usize]] { match self { ElementKind::Rod2 | ElementKind::Beam2 => &[&[0, 1]], + ElementKind::Rod3 => &[&[0, 1, 2]], ElementKind::Tri3 => &[&[0, 1], &[1, 2], &[2, 0]], ElementKind::Tri6 => &[&[0, 1, 3], &[1, 2, 4], &[2, 0, 5]], ElementKind::Quad4 => &[&[0, 1], &[1, 2], &[2, 3], &[3, 0]], @@ -357,7 +363,7 @@ mod tests { /// Reference coordinates with positive volume for each kind (mid-side nodes at edge midpoints). fn reference(kind: ElementKind) -> Vec<[f64; 3]> { let corners: Vec<[f64; 3]> = match kind { - ElementKind::Rod2 | ElementKind::Beam2 => vec![[0., 0., 0.], [1., 0., 0.]], + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => vec![[0., 0., 0.], [1., 0., 0.]], ElementKind::Tri3 | ElementKind::Tri6 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]], ElementKind::Quad4 | ElementKind::Quad8 => vec![[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.]], ElementKind::Tet4 | ElementKind::Tet10 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]], diff --git a/crates/model/src/sets.rs b/crates/model/src/sets.rs index 4ee451f..2df63c7 100644 --- a/crates/model/src/sets.rs +++ b/crates/model/src/sets.rs @@ -122,6 +122,7 @@ impl Model { let mut users = Vec::new(); users.extend(self.sections.iter().filter(|(_, s)| s.set == set).map(|(n, _)| format!("section `{n}`"))); users.extend(self.couplings.iter().filter(|(_, c)| c.set == set).map(|(n, _)| format!("coupling `{n}`"))); + users.extend(self.discrete_users(|t| t.set() == Some(set))); users .extend(self.boundary_conditions.iter().filter(|(_, b)| b.target().set() == Some(set)).map(|(n, _)| format!("boundary condition `{n}`"))); users.extend(self.loads.iter().filter(|(_, l)| l.set() == Some(set)).map(|(n, _)| format!("load `{n}`"))); diff --git a/crates/render/src/display.rs b/crates/render/src/display.rs index a66e294..d3630d0 100644 --- a/crates/render/src/display.rs +++ b/crates/render/src/display.rs @@ -79,7 +79,7 @@ pub struct DisplayMesh { /// counter-clockwise seen from outside for a positively oriented element. fn faces(kind: ElementKind) -> &'static [&'static [usize]] { match kind { - ElementKind::Rod2 | ElementKind::Beam2 => &[], + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => &[], ElementKind::Tri3 | ElementKind::Tri6 => &[&[0, 1, 2]], ElementKind::Quad4 | ElementKind::Quad8 => &[&[0, 1, 2, 3]], ElementKind::Tet4 | ElementKind::Tet10 => &[&[0, 2, 1], &[0, 1, 3], &[1, 2, 3], &[2, 0, 3]], diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 5f20529..91e7f46 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -12,7 +12,7 @@ use crate::constraints::{Budget, Reduction, coupling_rows}; use crate::eigen::{Operator, Wanted, lanczos}; use crate::ordering::nested_dissection; use crate::sparse::{Csr, NodeGraph, add_block, add_diagonal, add_scalar}; -use crate::system::{Mesh, inelastic_warnings, loads, prescribed}; +use crate::system::{Mesh, inelastic_warnings, loads, prescribed, translation_block}; use crate::{Error, Limits, Result}; /// Size and effort of a solution. @@ -147,6 +147,24 @@ fn assemble(mesh: &Mesh, graph: &NodeGraph, mass: Option) -> Result add_scalar(mm, layout, &t.nodes, kinds, &truss.mass(rho, t.area).map_err(at)?.data)?, + MassKind::Lumped => add_diagonal(mm, layout, &t.nodes, kinds, &truss.lumped_mass(rho, t.area).map_err(at)?)?, + } + } + } + mesh.discrete.add_stiffness(&mut k)?; + if let Some(mm) = m.as_mut() { + mesh.discrete.add_mass(mm)?; + } Ok(Assembled { k, m }) } @@ -184,6 +202,32 @@ fn split(full: &[f64], mesh: &Mesh) -> (Vec<[f64; 3]>, Vec<[f64; 6]>) { ((0..l.nodes()).map(|n| directions(full, l, n)).collect(), (0..mesh.points.len()).map(|k| directions(full, l, l.point_slot(k))).collect()) } +/// The matrix pattern of the mesh (node blocks of its elements, trusses and discrete elements) +/// and the order of its nodes for the factorization. +fn pattern(mesh: &Mesh) -> (NodeGraph, Vec) { + let graph = NodeGraph::new(mesh.block_count(), mesh.blocks()); + // The order covers the mesh nodes; reference point equations follow them. + let n = mesh.labels.len(); + let order = if graph.adjacency.len() == n { + nested_dissection(&mesh.xyz, &graph.adjacency) + } else { + let mesh_part: Vec> = graph.adjacency.iter().take(n).map(|row| row.iter().copied().filter(|b| *b < n).collect()).collect(); + nested_dissection(&mesh.xyz, &mesh_part) + }; + (graph, order) +} + +/// Warnings for what a step without velocities ignores. +fn ignored_dashpots(mesh: &Mesh, warnings: &mut Vec) { + let names = mesh.discrete.dashpots(); + if !names.is_empty() { + warnings.push(format!( + "dashpots {} act in dynamic steps only; this step ignores them", + names.iter().map(|n| format!("`{n}`")).collect::>().join(", ") + )); + } +} + /// Constraint rows, prescribed values and their elimination for a step. /// The model, its resolved view, the mesh and the step being solved. struct Setup<'a> { @@ -210,8 +254,8 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).collect(); warnings.extend(inelastic_warnings(model)); - let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice())); - let order = nested_dissection(&mesh.xyz, &graph.adjacency); + let (graph, order) = pattern(&mesh); + ignored_dashpots(&mesh, &mut warnings); let setup = Setup { model, resolved: &resolved, mesh: &mesh, step }; let red = reduction(&setup, true, &order, budget(limits, &mesh, &graph, 1), &mut warnings)?; let f = loads(model, &resolved, &mesh, step, &red.active)?; @@ -276,10 +320,28 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul } elements.push(ElementStress { label: e.label, kind: e.kind, points, nodal }); } + for t in &mesh.trusses { + let truss = t.truss()?; + let at = |err| Error::element(t.label, err); + let ue: Vec<[f64; 3]> = t.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect(); + let points = truss.stresses(t.young, &ue).map_err(at)?; + let nodal = truss.extrapolate(&points.iter().map(|p| p.stress).collect::>()).map_err(at)?; + for (n, s) in t.nodes.iter().zip(&nodal) { + if let (Some(acc), Some(c)) = (sum.get_mut(*n), count.get_mut(*n)) { + for k in 0..6 { + acc[k] += s[k]; + } + *c += 1; + } + } + elements.push(ElementStress { label: t.label, kind: t.kind, points, nodal }); + } + elements.sort_by_key(|e| e.label); let stress: Vec<[f64; 6]> = sum.iter().zip(&count).map(|(s, c)| s.map(|v| if *c > 0 { v / f64::from(*c) } else { 0.0 })).collect(); let von_mises = stress.iter().map(von_mises).collect(); let scale = applied_total.iter().chain(&reaction_total).fold(0.0f64, |m, v| m.max(v.abs())); - let imbalance = (0..3).map(|c| (applied_total[c] + reaction_total[c]).abs()).fold(0.0, f64::max); + let ground = mesh.discrete.ground_force(&u, &mesh); + let imbalance = (0..3).map(|c| (applied_total[c] + reaction_total[c] + ground[c]).abs()).fold(0.0, f64::max); if scale > 0.0 && imbalance > 1e-6 * scale { warnings.push(format!("loads and reactions differ by {imbalance:.3e} (relative {:.1e})", imbalance / scale)); } @@ -338,14 +400,23 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).collect(); warnings.extend(inelastic_warnings(model)); - let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice())); - let order = nested_dissection(&mesh.xyz, &graph.adjacency); + let (graph, order) = pattern(&mesh); + ignored_dashpots(&mesh, &mut warnings); let setup = Setup { model, resolved: &resolved, mesh: &mesh, step }; let red = reduction(&setup, false, &order, budget(limits, &mesh, &graph, 2), &mut warnings)?; let Assembled { k, m } = assemble(&mesh, &graph, Some(options.mass))?; let m = m.ok_or_else(|| Error::Internal("mass matrix".into()))?; // Every element's mass entries sum to its mass, once per direction. - let total_mass = m.val.iter().sum::() / if mesh.theory.is_some() { 2.0 } else { 3.0 }; + // (the mean over the directions for a point mass given per direction; rotary inertia is left + // out), added in storage order. + let translation = |d: usize| mesh.kind(d).and_then(DofKind::direction).is_some_and(|c| c < 3); + let mut translational = 0.0; + for d in (0..m.n).filter(|d| translation(*d)) { + for v in m.row(d).1 { + translational += v; + } + } + let total_mass = translational / mesh.translation_kinds().len() as f64; let neq = red.neq(); if neq == 0 { return Err(Error::Invalid("every degree of freedom is held; there is nothing to vibrate".into())); diff --git a/crates/solver/src/constraints.rs b/crates/solver/src/constraints.rs index fd5ee5a..202a295 100644 --- a/crates/solver/src/constraints.rs +++ b/crates/solver/src/constraints.rs @@ -191,10 +191,15 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) -> } /// The degrees of freedom that exist: every mesh node's, and the reference point directions -/// some coupling row involves. +/// some coupling row or discrete element involves. fn active_dofs(mesh: &Mesh, rows: &[Row]) -> Vec { let mut active = vec![false; mesh.dofs()]; active.iter_mut().take(mesh.mesh_dofs()).for_each(|a| *a = true); + for dof in mesh.discrete.dofs() { + if let Some(a) = active.get_mut(dof) { + *a = true; + } + } for row in rows { for (dof, c) in &row.terms { if *c != 0.0 diff --git a/crates/solver/src/discrete.rs b/crates/solver/src/discrete.rs new file mode 100644 index 0000000..fe21b42 --- /dev/null +++ b/crates/solver/src/discrete.rs @@ -0,0 +1,198 @@ +//! Point masses, rotary inertia, springs and dashpots in the implicit system. +//! +//! The model resolves each definition to its nodes ([`modelcraft_model::Resolved::masses`] and +//! `springs`); here each copy is turned into global degrees of freedom and its matrix: +//! +//! - a point mass adds m_x, m_y, m_z on the node's translations and its inertia tensor J (3 × 3, +//! off-diagonal products included) on its rotations, the same in the consistent and the lumped +//! mass matrix; +//! - a spring with stiffness k and unit direction coefficients c on each end's six degrees of +//! freedom adds k c cᵀ to each end and −k c cᵀ between them (Cook, Malkus, Plesha, Witt 2002, +//! §2.2): the elongation is cᵀ(u_second − u_first); +//! - a dashpot has the same form on velocities; static and frequency steps have no velocities, so +//! they ignore dashpots with a warning (damping enters dynamic steps). +//! +//! What each copy needs at its nodes ([`needs`]) goes into the mesh's [`modelcraft_fem::DofLayout`] +//! like an element's declared kinds; [`dof_of`] then finds its unknowns there. + +use modelcraft_fem::{DofKind, DofSet}; +use modelcraft_model::{Resolved, ResolvedSpring}; + +use crate::sparse::{Csr, add_entry}; +use crate::system::Mesh; +use crate::{Error, Result}; + +/// A point mass on global degrees of freedom. +#[derive(Debug, Clone, PartialEq)] +pub struct MassAt { + pub name: String, + /// Translations along X, Y, Z: (degree of freedom, mass). + pub translations: Vec<(usize, f64)>, + /// The rotations' degrees of freedom and the inertia tensor on them (empty without inertia). + pub rotations: Vec, + pub inertia: [[f64; 3]; 3], +} + +/// A spring or dashpot on global degrees of freedom: the elongation is Σ c u over `terms`, whose +/// coefficients are the direction's on the second end and their negatives on the first. +#[derive(Debug, Clone, PartialEq)] +pub struct SpringAt { + pub name: String, + pub dashpot: bool, + /// Connected to the ground (one end only). + pub grounded: bool, + pub terms: Vec<(usize, f64)>, + pub value: f64, +} + +/// Every discrete element of the model on the mesh's degrees of freedom. +#[derive(Debug, Clone, Default, PartialEq)] +pub struct Discrete { + pub masses: Vec, + pub springs: Vec, +} + +/// The kinds each discrete copy needs at the labels it acts on: translations for a mass, +/// rotations for an inertia, translations or rotations for a spring or dashpot. +/// `translations` are the model's (X and Y in a 2D model). +pub fn needs(resolved: &Resolved, translations: DofSet) -> impl Iterator + '_ { + let masses = resolved.masses.iter().flat_map(move |m| { + let mass = m.mass.iter().any(|c| *c != 0.0).then_some((m.node, translations)); + let inertia = m.inertia.iter().any(|c| *c != 0.0).then_some((m.node, DofSet::ROTATIONS)); + mass.into_iter().chain(inertia) + }); + // A node of a 3D model moves in all three directions (a spring along X does not take Y and + // Z away), so a spring asks for every translation, or every rotation. + let springs = resolved.springs.iter().flat_map(move |s| { + let kinds = if s.direction.is_rotational() { DofSet::ROTATIONS } else { translations }; + std::iter::once((s.first, kinds)).chain(s.second.map(|b| (b, kinds))) + }); + masses.chain(springs) +} + +/// The global degree of freedom of direction `d` (0–2 translations, 3–5 rotations) of a node or +/// reference point label, or `None` when the node has no such direction. +pub fn dof_of(mesh: &Mesh, label: u32, d: usize) -> Option { + mesh.layout.dof(mesh.layout.slot(label)?, DofKind::from_direction(d)?) +} + +fn missing(what: &str, name: &str, label: u32, d: usize) -> Error { + let direction = DofKind::from_direction(d).map_or("?", DofKind::describe); + Error::Invalid(format!("{what} `{name}`: node {label} has no degree of freedom in {direction}")) +} + +/// The spring's terms at one end, with the direction's coefficients times `sign`. +fn ends(mesh: &Mesh, s: &ResolvedSpring, label: u32, sign: f64) -> Result> { + s.direction + .coefficients() + .iter() + .enumerate() + .filter(|(_, c)| **c != 0.0) + .map(|(d, c)| dof_of(mesh, label, d).map(|dof| (dof, sign * c)).ok_or_else(|| missing(s.kind(), &s.name, label, d))) + .collect() +} + +impl Discrete { + /// The resolved discrete elements on the mesh's degrees of freedom. + pub fn new(mesh: &Mesh, resolved: &Resolved) -> Result { + let mut out = Discrete::default(); + for m in &resolved.masses { + // The model's translations only (a 2D model's masses act in its plane). + let planar = mesh.translation_kinds(); + let translations: Vec<(usize, f64)> = (0..3) + .filter(|d| m.mass[*d] != 0.0 && DofKind::translation(*d).is_some_and(|k| planar.contains(k))) + .map(|d| dof_of(mesh, m.node, d).map(|dof| (dof, m.mass[d])).ok_or_else(|| missing("point mass", &m.name, m.node, d))) + .collect::>()?; + let inertia = modelcraft_model::discrete::inertia_matrix(m.inertia); + let rotations = if m.inertia.iter().any(|c| *c != 0.0) { + (3..6).map(|d| dof_of(mesh, m.node, d).ok_or_else(|| missing("point mass", &m.name, m.node, d))).collect::>()? + } else { + Vec::new() + }; + out.masses.push(MassAt { name: m.name.clone(), translations, rotations, inertia }); + } + for s in &resolved.springs { + let mut terms = ends(mesh, s, s.first, -1.0)?; + if let Some(b) = s.second { + terms.extend(ends(mesh, s, b, 1.0)?); + } + out.springs.push(SpringAt { name: s.name.clone(), dashpot: s.dashpot, grounded: s.second.is_none(), terms, value: s.value }); + } + Ok(out) + } + + /// The groups of degrees of freedom each element couples, for the matrix pattern. + pub fn couplings(&self) -> impl Iterator> + '_ { + let masses = self.masses.iter().map(|m| m.translations.iter().map(|t| t.0).chain(m.rotations.iter().copied()).collect()); + masses.chain(self.springs.iter().map(|s| s.terms.iter().map(|t| t.0).collect())) + } + + /// Every degree of freedom a discrete element acts on. + pub fn dofs(&self) -> impl Iterator + '_ { + self.couplings().flatten() + } + + /// Names of the dashpots, which static and frequency steps ignore. + pub fn dashpots(&self) -> Vec<&str> { + let mut names: Vec<&str> = self.springs.iter().filter(|s| s.dashpot).map(|s| s.name.as_str()).collect(); + names.dedup(); + names + } + + /// Add the springs' stiffness. + pub fn add_stiffness(&self, k: &mut Csr) -> Result<()> { + for s in self.springs.iter().filter(|s| !s.dashpot) { + for (i, ci) in &s.terms { + for (j, cj) in &s.terms { + add_entry(k, *i, *j, s.value * ci * cj)?; + } + } + } + Ok(()) + } + + /// Add the point masses and rotary inertia. + pub fn add_mass(&self, m: &mut Csr) -> Result<()> { + for p in &self.masses { + for (d, v) in &p.translations { + add_entry(m, *d, *d, *v)?; + } + for (a, da) in p.rotations.iter().enumerate() { + for (b, db) in p.rotations.iter().enumerate() { + let v = p.inertia[a][b]; + if v != 0.0 { + add_entry(m, *da, *db, v)?; + } + } + } + } + Ok(()) + } + + /// The sum of the forces springs to the ground exert on the model, along X, Y and Z (the + /// loads they carry to the ground, besides the boundary conditions' reactions). + pub fn ground_force(&self, u: &[f64], mesh: &Mesh) -> [f64; 3] { + let mut total = [0.0; 3]; + for s in self.springs.iter().filter(|s| !s.dashpot && s.grounded) { + // The spring's force on the model: −k (cᵀu) c over its terms. + let stretch: f64 = s.terms.iter().map(|(d, c)| c * u.get(*d).copied().unwrap_or(0.0)).sum(); + for (d, c) in &s.terms { + if let Some(t) = mesh.kind(*d).and_then(DofKind::direction).filter(|d| *d < 3).and_then(|d| total.get_mut(d)) { + *t -= s.value * stretch * c; + } + } + } + total + } + + /// Gravity on the point masses: m_d g_d on each translation. + pub fn add_gravity(&self, f: &mut [f64], g: [f64; 3], mesh: &Mesh) { + for p in &self.masses { + for (dof, m) in &p.translations { + if let (Some(x), Some(a)) = (f.get_mut(*dof), mesh.kind(*dof).and_then(DofKind::direction).and_then(|d| g.get(d))) { + *x += m * a; + } + } + } + } +} diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index 56fffc0..50db5b9 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -25,6 +25,7 @@ pub mod analysis; pub mod cholesky; pub mod constraints; +pub mod discrete; pub mod eigen; pub mod ordering; pub mod sparse; diff --git a/crates/solver/src/sparse.rs b/crates/solver/src/sparse.rs index c078d32..e57b5a3 100644 --- a/crates/solver/src/sparse.rs +++ b/crates/solver/src/sparse.rs @@ -103,10 +103,10 @@ impl NodeGraph { self.adjacency.iter().enumerate().map(|(a, row)| layout.range(a).len() * row.iter().map(|b| layout.range(*b).len()).sum::()).sum() } - /// A zero matrix over the mesh nodes' degrees of freedom (layout slots `0..nodes`) on this - /// pattern. + /// A zero matrix over the degrees of freedom of the graph's layout slots (the mesh nodes + /// `0..nodes`, and reference points when the graph has rows for them) on this pattern. pub fn zero_matrix(&self, layout: &DofLayout) -> Csr { - let n = layout.node_dofs(); + let n = self.adjacency.len().checked_sub(1).map_or(0, |last| layout.range(last).end); let mut ptr = Vec::with_capacity(n + 1); ptr.push(0); let mut idx = Vec::with_capacity(self.entries(layout)); @@ -210,6 +210,14 @@ pub fn add_diagonal(k: &mut Csr, layout: &DofLayout, nodes: &[usize], kinds: Dof Ok(()) } +/// Add `v` to entry (`i`, `j`), which must be in the pattern. +pub fn add_entry(k: &mut Csr, i: usize, j: usize, v: f64) -> Result<()> { + let (start, end) = (k.ptr.get(i).copied().unwrap_or(0), k.ptr.get(i + 1).copied().unwrap_or(0)); + let at = k.idx.get(start..end).and_then(|cols| cols.binary_search(&j).ok()).ok_or_else(|| Error::Internal("entry outside the pattern".into()))?; + add_at(k, start + at, v); + Ok(()) +} + #[cfg(test)] mod tests { use super::*; diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 1179394..943eaab 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -1,5 +1,6 @@ -//! From the resolved model to the discrete system: the solid elements with their materials and -//! formulations, the coupled reference points, the degrees of freedom, prescribed values and the +//! From the resolved model to the discrete system: the solid, 2D and truss elements with their +//! materials and formulations, point masses, springs and dashpots ([`crate::discrete`]), the +//! reference points couplings or discrete elements use, the degrees of freedom, prescribed values and the //! load vector of a step. //! //! Degrees of freedom ([`DofLayout`]): the kinds each mesh node's elements declare (three @@ -8,6 +9,7 @@ use modelcraft_fem::element::{Formulation, Solid}; use modelcraft_fem::plane::{Plane, Theory}; +use modelcraft_fem::truss::{Truss, young_modulus}; use modelcraft_fem::{DofKind, DofLayout, DofSet, Elasticity, ElementDofs, PointStress, elasticity, face, plane}; use modelcraft_model::nodal::edge_pressure; use std::collections::BTreeSet; @@ -117,6 +119,42 @@ impl Elem { } } +/// A truss element (rod2, rod3) ready for the kernel. +#[derive(Debug, Clone)] +pub struct Bar { + pub label: u32, + pub kind: ElementKind, + /// Indices into [`Mesh::labels`] (layout slots). + pub nodes: Vec, + pub xyz: Vec<[f64; 3]>, + pub young: f64, + pub area: f64, + pub density: Option, +} + +impl Bar { + pub fn truss(&self) -> Result> { + Truss::new(self.kind, &self.xyz).map_err(|e| Error::element(self.label, e)) + } +} + +/// The translations of a model's nodes: X and Y in a 2D model, else X, Y and Z. Trusses and +/// translational discrete elements declare these. +pub fn translation_kinds(planar: bool) -> DofSet { + if planar { DofSet::PLANE } else { DofSet::TRANSLATIONS } +} + +/// Keep the rows and columns of a matrix over three translations per node (node by node) that +/// `kinds` has: the in-plane block of a truss in a 2D model, or the matrix itself. +pub fn translation_block(data: &[f64], nodes: usize, kinds: DofSet) -> Vec { + if kinds == DofSet::TRANSLATIONS { + return data.to_vec(); + } + let keep: Vec = (0..3 * nodes).filter(|i| DofKind::translation(i % 3).is_some_and(|k| kinds.contains(k))).collect(); + let width = 3 * nodes; + keep.iter().flat_map(|i| keep.iter().map(move |j| data.get(i * width + j).copied().unwrap_or(0.0))).collect() +} + /// A reference point some coupling ties to the mesh. #[derive(Debug, Clone, PartialEq)] pub struct Point { @@ -138,6 +176,9 @@ pub struct Mesh { pub layout: DofLayout, /// The 2D theory when this is a 2D model (a model is 2D as a whole). pub theory: Option, + pub trusses: Vec, + /// Point masses, springs and dashpots ([`crate::discrete`]). + pub discrete: crate::discrete::Discrete, } impl Mesh { @@ -145,6 +186,15 @@ impl Mesh { self.labels.binary_search(&label).ok() } + pub fn truss(&self, label: u32) -> Option<&Bar> { + self.trusses.binary_search_by_key(&label, |e| e.label).ok().and_then(|i| self.trusses.get(i)) + } + + /// The kinds trusses and translational discrete elements carry at their nodes. + pub fn translation_kinds(&self) -> DofSet { + translation_kinds(self.theory.is_some()) + } + pub fn element(&self, label: u32) -> Option<&Elem> { self.elements.binary_search_by_key(&label, |e| e.label).ok().and_then(|i| self.elements.get(i)) } @@ -188,17 +238,32 @@ impl Mesh { /// Build it from the model and its resolved view, refusing what the solver cannot represent. pub fn new(model: &Model, resolved: &Resolved, max_dofs: usize) -> Result { let mut elements = Vec::new(); + let mut trusses = Vec::new(); // A model is 2D as a whole (the model check refuses a mix). let planar = model.sections().values().find_map(|s| Theory::from_section(&s.kind)); + if matches!(planar, Some(Theory::Axisymmetric)) + && let Some(what) = model + .sections() + .values() + .any(|s| matches!(s.kind, SectionKind::Rod { .. })) + .then_some("trusses") + .or_else(|| (!model.masses().is_empty() || !model.springs().is_empty()).then_some("point masses, springs and dashpots")) + { + return Err(Error::NotSupported(format!("{what} in axisymmetric models"))); + } for (name, section) in model.sections() { let theory = Theory::from_section(§ion.kind); - if !matches!(section.kind, SectionKind::Solid {}) && theory.is_none() { + let rod = match section.kind { + SectionKind::Rod { area } => Some(area), + _ => None, + }; + if !matches!(section.kind, SectionKind::Solid {}) && theory.is_none() && rod.is_none() { return Err(Error::NotSupported(format!( - "{} sections (section `{name}`); the solver has solid and 2D (plane stress, plane strain, axisymmetric) elements so far", + "{} sections (section `{name}`); the solver has solid, 2D (plane stress, plane strain, axisymmetric) and truss elements so far", section.kind.name() ))); } - if theory.is_some() != planar.is_some() { + if theory.is_some() != planar.is_some() && rod.is_none() { return Err(Error::NotSupported("2D and solid sections in one model".into())); } let material = model @@ -206,8 +271,17 @@ impl Mesh { .ok_or_else(|| Error::Invalid(format!("section `{name}`: material `{}` is missing", section.material)))?; let elastic = material.elastic.as_ref().ok_or_else(|| Error::Invalid(format!("material `{}` has no elastic constants", section.material)))?; - let d = elasticity(elastic).map_err(|e| Error::Invalid(format!("material `{}`: {e}", section.material)))?; let set = resolved.sets.get(§ion.set).ok_or_else(|| Error::Invalid(format!("set `{}` is not resolved", section.set)))?; + if let Some(area) = rod { + let young = + young_modulus(elastic).map_err(|e| Error::Invalid(format!("material `{}` of section `{name}`: {e}", section.material)))?; + for &label in &set.elements { + let el = model.element(label).ok_or_else(|| Error::Invalid(format!("element {label} is missing")))?; + trusses.push((label, el.kind, el.nodes.clone(), young, area, material.density)); + } + continue; + } + let d = elasticity(elastic).map_err(|e| Error::Invalid(format!("material `{}`: {e}", section.material)))?; for &label in &set.elements { let el = model.element(label).ok_or_else(|| Error::Invalid(format!("element {label} is missing")))?; let (formulation, dofs) = match theory { @@ -222,10 +296,19 @@ impl Mesh { return Err(Error::NotSupported(format!("couplings in 2D models (coupling `{name}`)"))); } elements.sort_by_key(|e| e.0); + trusses.sort_by_key(|e| e.0); + let translations = translation_kinds(planar.is_some()); + // What discrete elements need: at mesh nodes (some only they use), and at reference + // points, which then get their six degrees of freedom like coupled ones. + let needs: Vec<(u32, DofSet)> = crate::discrete::needs(resolved, translations).collect(); let mut labels: Vec = elements.iter().flat_map(|e| e.2.iter().copied()).collect(); + labels.extend(trusses.iter().flat_map(|e| e.2.iter().copied())); + labels.extend(needs.iter().map(|n| n.0).filter(|l| model.node(*l).is_some())); labels.sort_unstable(); labels.dedup(); - let coupled: BTreeSet<&str> = model.couplings().values().map(|c| c.reference_point.as_str()).collect(); + let mut coupled: BTreeSet<&str> = model.couplings().values().map(|c| c.reference_point.as_str()).collect(); + let discrete_labels: BTreeSet = needs.iter().map(|n| n.0).collect(); + coupled.extend(model.reference_points().iter().filter(|(_, p)| discrete_labels.contains(&p.label)).map(|(n, _)| n.as_str())); let points: Vec = coupled .into_iter() .map(|name| { @@ -237,6 +320,12 @@ impl Mesh { for e in &elements { layout.element(&e.2, e.7); } + for e in &trusses { + layout.element(&e.2, modelcraft_fem::ElementDofs::new(translations)); + } + for (label, kinds) in needs.iter().filter(|n| model.node(n.0).is_some()) { + layout.node(*label, *kinds); + } for p in &points { layout.point(p.label, DofSet::STRUCTURAL); } @@ -249,14 +338,46 @@ impl Mesh { .iter() .map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| Error::Invalid(format!("node {l} is missing")))) .collect::>()?; - let mut mesh = Mesh { labels, xyz, elements: Vec::with_capacity(elements.len()), points, layout, theory: planar }; + let mut mesh = Mesh { + labels, + xyz, + elements: Vec::with_capacity(elements.len()), + points, + layout, + theory: planar, + trusses: Vec::with_capacity(trusses.len()), + discrete: crate::discrete::Discrete::default(), + }; for (label, kind, nodes, formulation, d, density, theory, _) in elements { let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); mesh.elements.push(Elem { label, kind, nodes: idx, xyz: exyz, formulation, d, density, theory }); } + for (label, kind, nodes, young, area, density) in trusses { + let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); + let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); + mesh.trusses.push(Bar { label, kind, nodes: idx, xyz: exyz, young, area, density }); + } + mesh.discrete = crate::discrete::Discrete::new(&mesh, resolved)?; Ok(mesh) } + + /// The layout slots (mesh nodes, reference points) every element and discrete element + /// couples, for the matrix pattern. + pub fn blocks(&self) -> impl Iterator> { + let elements = self.elements.iter().map(|e| e.nodes.as_slice()).chain(self.trusses.iter().map(|t| t.nodes.as_slice())); + elements.map(std::borrow::Cow::Borrowed).chain(self.discrete.couplings().map(|dofs| { + let mut slots: Vec = dofs.iter().filter_map(|d| self.layout.locate(*d).map(|(slot, _)| slot)).collect(); + slots.dedup(); + std::borrow::Cow::Owned(slots) + })) + } + + /// The slots of the matrix pattern: the mesh nodes, and the reference points too when a + /// discrete element acts on them (couplings tie points through constraints, not the matrix). + pub fn block_count(&self) -> usize { + if self.discrete.couplings().flatten().any(|d| d >= self.mesh_dofs()) { self.layout.slots() } else { self.labels.len() } + } } /// The prescribed value of every degree of freedom (`None` = free) under a step's boundary @@ -279,8 +400,9 @@ pub fn prescribed( if let Target::ReferencePoint(rp) = bc.target() && mesh.point(rp).is_none() { - warnings - .push(format!("boundary condition `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh; it holds nothing")); + warnings.push(format!( + "boundary condition `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh and no spring or point mass uses; it holds nothing" + )); continue; } for slot in resolved.target_nodes(bc.target()).iter().filter_map(|l| mesh.layout.slot(*l)) { @@ -365,7 +487,9 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ if let Some(rp) = load.reference_point() && mesh.point(rp).is_none() { - return Err(Error::Invalid(format!("load `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh"))); + return Err(Error::Invalid(format!( + "load `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh and no spring or point mass uses" + ))); } // Moments act on the rotations, which only reference points have. let rotation = matches!(load, Load::Moment { .. }); @@ -433,15 +557,36 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ let rho = e.density.ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} has no density", e.label)))?; body(&mut f, mesh, e, acceleration.map(|a| rho * a))?; } + for t in &mesh.trusses { + if chosen.is_some_and(|c| c.binary_search(&t.label).is_err()) { + continue; + } + let rho = t.density.ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} has no density", t.label)))?; + bar_body(&mut f, mesh, t, acceleration.map(|a| rho * a))?; + } + // Point masses take gravity on the whole model. + if set.is_none() { + mesh.discrete.add_gravity(&mut f, *acceleration, mesh); + } } Load::BodyForce { set, vector } => { for label in resolved.sets.get(set).map(|r| r.elements.as_slice()).unwrap_or_default() { + if let Some(t) = mesh.truss(*label) { + bar_body(&mut f, mesh, t, *vector)?; + continue; + } let e = mesh.element(*label).ok_or_else(|| Error::Invalid(format!("load `{name}`: element {label} is not a solid element")))?; body(&mut f, mesh, e, *vector)?; } } Load::Centrifugal { angular_velocity, point, axis, set } => { let chosen: Option<&[u32]> = set.as_ref().map(|s| resolved.sets.get(s).map(|r| r.elements.as_slice()).unwrap_or_default()); + if let Some(t) = mesh.trusses.iter().find(|t| chosen.is_none_or(|c| c.binary_search(&t.label).is_ok())) { + return Err(Error::NotSupported(format!("centrifugal loads on trusses (load `{name}`, element {})", t.label))); + } + if set.is_none() && !mesh.discrete.masses.is_empty() { + return Err(Error::NotSupported(format!("centrifugal loads on point masses (load `{name}`)"))); + } for e in &mesh.elements { if chosen.is_some_and(|c| c.binary_search(&e.label).is_err()) { continue; @@ -461,6 +606,15 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ Ok(f) } +/// A force `b` per volume on a truss. +fn bar_body(f: &mut [f64], mesh: &Mesh, t: &Bar, b: [f64; 3]) -> Result<()> { + let forces = t.truss()?.body_force(b, t.area).map_err(|err| Error::element(t.label, err))?; + for (node, v) in t.nodes.iter().zip(forces) { + add(f, &translations(mesh, *node)?, v); + } + Ok(()) +} + fn body(f: &mut [f64], mesh: &Mesh, e: &Elem, b: [f64; 3]) -> Result<()> { let forces = e.kernel()?.body_force(b).map_err(|err| Error::element(e.label, err))?; for (node, v) in e.nodes.iter().zip(forces) { diff --git a/crates/solver/tests/one_d.rs b/crates/solver/tests/one_d.rs new file mode 100644 index 0000000..ac099b0 --- /dev/null +++ b/crates/solver/tests/one_d.rs @@ -0,0 +1,652 @@ +//! Trusses, point masses, rotary inertia and springs against closed-form solutions. Run with +//! `--nocapture` to see the errors. +#![allow(clippy::needless_range_loop)] +//! +//! - Statically determinate trusses (a plane Pratt truss and a space tripod): the member forces +//! are those of the method of joints, whatever the members' stiffness. +//! - The three-bar truss (statically indeterminate): N = P cos²θ / (1 + 2 cos³θ) in the inclined +//! bars, P / (1 + 2 cos³θ) in the vertical one. +//! - rod3: constant-stress patch with off-centre middle nodes; a bar hanging under its own +//! weight (quadratic displacement, exact at the nodes); curved members are checked in +//! `modelcraft_fem::truss` (a ring's radial expansion gives the hoop strain δ/R exactly). +//! - Spring–mass ω = √(k/m), springs along three skew directions, the two-mass chain, a bar with +//! a tip mass ω = √(EA/(LM)), the torsional pendulum ω = √(k_t/J). + +use std::f64::consts::PI; + +use modelcraft_fem::truss::axial; +use modelcraft_model::*; +use modelcraft_solver::{FrequencyOptions, FrequencyResult, Limits, MassKind, StaticResult, solve_frequency, solve_static}; + +const ZERO_Z: [Option; 6] = [None, None, Some(0.0), None, None, None]; + +/// A model built point by point. +struct Build { + m: Model, + step: &'static str, +} + +impl Build { + fn new(procedure: Procedure) -> Build { + let mut m = Model::default(); + m.add_step(Step::new("Step", procedure), None).unwrap(); + Build { m, step: "Step" } + } + + fn statics() -> Build { + Build::new(Procedure::Static {}) + } + + fn modes(n: u32) -> Build { + Build::new(Procedure::Frequency { modes: n, min_hz: None, max_hz: None }) + } + + fn node(&mut self, label: u32, xyz: [f64; 3]) { + self.m.set_node(label, xyz).unwrap(); + } + + fn nodes(&mut self, name: &str, nodes: &[u32]) -> Target { + self.m.add_set(name, Set { nodes: nodes.to_vec(), ..Set::default() }).unwrap(); + Target::Set(name.into()) + } + + /// Members with one area and material. + fn members(&mut self, name: &str, kind: ElementKind, members: &[(u32, Vec)], area: f64, young: f64, density: Option) { + for (label, nodes) in members { + self.m.set_element(*label, kind, nodes.clone()).unwrap(); + } + self.m.add_set(name, Set { elements: members.iter().map(|m| m.0).collect(), ..Set::default() }).unwrap(); + self.m.add_material(name, Material { elastic: Some(Elastic::Isotropic { young, poisson: 0.3 }), density, ..Material::default() }).unwrap(); + self.m + .add_section( + name, + Section { set: name.into(), material: name.into(), kind: SectionKind::Rod { area }, formulation: Formulation::default() }, + ) + .unwrap(); + } + + fn bc(&mut self, name: &str, target: Target, values: [Option; 6]) { + self.m.add_boundary_condition(name, BoundaryCondition::Displacement { target, values }).unwrap(); + self.m.set_boundary_condition_active(self.step, name, true).unwrap(); + } + + fn fix(&mut self, name: &str, target: Target) { + self.m.add_boundary_condition(name, BoundaryCondition::Fixed { target }).unwrap(); + self.m.set_boundary_condition_active(self.step, name, true).unwrap(); + } + + fn load(&mut self, name: &str, load: Load) { + self.m.add_load(name, load).unwrap(); + self.m.set_load_active(self.step, name, true).unwrap(); + } + + fn force(&mut self, name: &str, node: u32, vector: [f64; 3]) { + let target = self.nodes(name, &[node]); + self.load(name, Load::Force { target, vector, distribution: Distribution::Total }); + } + + fn spring(&mut self, name: &str, first: Target, second: Option, direction: Direction, stiffness: f64) { + self.m.add_spring(name, Spring::Spring { first, second, direction, stiffness }).unwrap(); + } + + fn mass(&mut self, name: &str, target: Target, mass: Option, inertia: Option<[f64; 6]>) { + self.m.add_mass(name, PointMass { target, mass, inertia }).unwrap(); + } + + fn solve(&self) -> StaticResult { + solve_static(&self.m, None, &Limits::default()).unwrap() + } + + fn modes_of(&self, mass: MassKind) -> FrequencyResult { + solve_frequency(&self.m, None, &FrequencyOptions { mass, shift: None }, &Limits::default()).unwrap() + } +} + +fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[0] - b[0], a[1] - b[1], a[2] - b[2]] +} + +fn unit(a: [f64; 3]) -> [f64; 3] { + let l = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt(); + a.map(|c| c / l) +} + +/// The axial force of each member at its first integration point: N = A tᵀ σ t. +fn member_forces(m: &Model, r: &StaticResult, area_of: impl Fn(u32) -> f64) -> Vec<(u32, f64)> { + r.elements + .iter() + .map(|e| { + let el = m.element(e.label).unwrap(); + let t = unit(sub(m.node(el.nodes[1]).unwrap().xyz, m.node(el.nodes[0]).unwrap().xyz)); + (e.label, area_of(e.label) * axial(&e.points[0].stress, t, false)) + }) + .collect() +} + +/// Gaussian elimination with partial pivoting (the method of joints' equations). +fn solve_dense(mut a: Vec>, mut b: Vec) -> Vec { + let n = b.len(); + for c in 0..n { + let p = (c..n).max_by(|i, j| a[*i][c].abs().total_cmp(&a[*j][c].abs())).unwrap(); + a.swap(c, p); + b.swap(c, p); + assert!(a[c][c].abs() > 1e-12, "the truss is not statically determinate"); + for r in c + 1..n { + let f = a[r][c] / a[c][c]; + for k in c..n { + a[r][k] -= f * a[c][k]; + } + b[r] -= f * b[c]; + } + } + let mut x = vec![0.0; n]; + for r in (0..n).rev() { + x[r] = (b[r] - (r + 1..n).map(|k| a[r][k] * x[k]).sum::()) / a[r][r]; + } + x +} + +/// Member forces by the method of joints: at every joint and direction, Σ N t + F + R = 0, with +/// the reactions `supports` (node, direction) as unknowns too. +fn method_of_joints( + xyz: &[(u32, [f64; 3])], + members: &[(u32, [u32; 2])], + supports: &[(u32, usize)], + loads: &[(u32, [f64; 3])], + dims: usize, +) -> Vec<(u32, f64)> { + let row = |node: u32, d: usize| xyz.iter().position(|(l, _)| *l == node).unwrap() * dims + d; + let pos = |node: u32| xyz.iter().find(|(l, _)| *l == node).unwrap().1; + let n = xyz.len() * dims; + assert_eq!(members.len() + supports.len(), n, "determinate"); + let mut a = vec![vec![0.0; n]; n]; + let mut b = vec![0.0; n]; + for (k, (_, [i, j])) in members.iter().enumerate() { + let t = unit(sub(pos(*j), pos(*i))); + for d in 0..dims { + // Tension pulls joint i towards j and j towards i. + a[row(*i, d)][k] += t[d]; + a[row(*j, d)][k] -= t[d]; + } + } + for (k, (node, d)) in supports.iter().enumerate() { + a[row(*node, *d)][members.len() + k] = 1.0; + } + for (node, f) in loads { + for d in 0..dims { + b[row(*node, d)] -= f[d]; + } + } + let x = solve_dense(a, b); + members.iter().zip(x).map(|((l, _), n)| (*l, n)).collect() +} + +fn worst_relative(fe: &[(u32, f64)], exact: &[(u32, f64)]) -> f64 { + let scale = exact.iter().fold(0.0f64, |m, (_, n)| m.max(n.abs())); + exact.iter().map(|(l, n)| (fe.iter().find(|(m, _)| m == l).unwrap().1 - n).abs() / scale).fold(0.0, f64::max) +} + +/// A plane Pratt truss of four panels, pinned at one end and on a roller at the other, loaded at +/// the bottom joints; chords and web members have different areas. The out-of-plane direction is +/// held. Member forces equal the method of joints'. +#[test] +fn pratt_truss_member_forces_by_the_method_of_joints() { + let (a, h, p) = (3.0, 4.0, 10.0); + let mut b = Build::statics(); + let mut xyz = Vec::new(); + for i in 0..5u32 { + xyz.push((1 + i, [a * f64::from(i), 0.0, 0.0])); + } + for i in 1..4u32 { + xyz.push((10 + i, [a * f64::from(i), h, 0.0])); + } + for (l, x) in &xyz { + b.node(*l, *x); + } + let chords: Vec<(u32, [u32; 2])> = vec![(1, [1, 2]), (2, [2, 3]), (3, [3, 4]), (4, [4, 5]), (5, [11, 12]), (6, [12, 13])]; + let web: Vec<(u32, [u32; 2])> = vec![(7, [2, 11]), (8, [3, 12]), (9, [4, 13]), (10, [1, 11]), (11, [13, 5]), (12, [11, 3]), (13, [13, 3])]; + let as_members = |list: &[(u32, [u32; 2])]| list.iter().map(|(l, n)| (*l, n.to_vec())).collect::>(); + b.members("chords", ElementKind::Rod2, &as_members(&chords), 2.0, 2.0e5, None); + b.members("web", ElementKind::Rod2, &as_members(&web), 0.7, 7.0e4, None); + let all: Vec = xyz.iter().map(|x| x.0).collect(); + let plane = b.nodes("all_nodes", &all); + b.bc("plane", plane, ZERO_Z); + let pin = b.nodes("pin", &[1]); + b.bc("pin", pin, [Some(0.0), Some(0.0), None, None, None, None]); + let roller = b.nodes("roller", &[5]); + b.bc("roller", roller, [None, Some(0.0), None, None, None, None]); + let loads = [(2, [0.0, -p, 0.0]), (3, [0.0, -1.5 * p, 0.0]), (4, [0.2 * p, -p, 0.0])]; + for (node, f) in loads { + b.force(&format!("P{node}"), node, f); + } + let r = b.solve(); + let fe = member_forces(&b.m, &r, |l| if l <= 6 { 2.0 } else { 0.7 }); + let members: Vec<(u32, [u32; 2])> = chords.iter().chain(&web).copied().collect(); + let exact = method_of_joints(&xyz, &members, &[(1, 0), (1, 1), (5, 1)], &loads, 2); + let err = worst_relative(&fe, &exact); + eprintln!("Pratt truss: worst member force error {err:.2e} of the largest"); + assert!(err < 1e-10, "{fe:?} vs {exact:?}"); + // Reactions balance the loads. + for d in 0..2 { + assert!((r.applied_total[d] + r.reaction_total[d]).abs() < 1e-9 * p); + } +} + +/// A space tripod: three members from an apex to three fixed points, a skew load at the apex. +#[test] +fn space_tripod_member_forces() { + let (rad, h) = (2.0, 3.0); + let mut b = Build::statics(); + let apex = (4, [0.1, -0.2, h]); + let feet: Vec<(u32, [f64; 3])> = (0..3) + .map(|i| (1 + i as u32, [rad * (PI / 2.0 + 2.0 * PI * i as f64 / 3.0).cos(), rad * (PI / 2.0 + 2.0 * PI * i as f64 / 3.0).sin(), 0.0])) + .collect(); + for (l, x) in feet.iter().chain([&apex]) { + b.node(*l, *x); + } + let members: Vec<(u32, Vec)> = (1..=3).map(|i| (i, vec![i, 4])).collect(); + b.members("legs", ElementKind::Rod2, &members, 1.5, 1.0e5, None); + let feet_set = b.nodes("feet", &[1, 2, 3]); + b.fix("feet", feet_set); + let load = [3.0, -2.0, -10.0]; + b.force("P", 4, load); + let r = b.solve(); + let fe = member_forces(&b.m, &r, |_| 1.5); + // Equilibrium of the apex: Σ N_i t_i + P = 0 with t_i from the apex to foot i (tension + // pulls the apex towards the foot). + let t: Vec<[f64; 3]> = feet.iter().map(|(_, x)| unit(sub(*x, apex.1))).collect(); + let a = (0..3).map(|d| (0..3).map(|i| t[i][d]).collect()).collect(); + let exact: Vec<(u32, f64)> = (1..=3).zip(solve_dense(a, load.map(|c| -c).to_vec())).collect(); + let err = worst_relative(&fe, &exact); + eprintln!("space tripod: worst member force error {err:.2e}"); + assert!(err < 1e-10, "{fe:?} vs {exact:?}"); +} + +/// The statically indeterminate three-bar truss: a vertical bar of length L and two at ±θ +/// meeting at a loaded joint. +#[test] +fn three_bar_truss() { + let (l, theta, p, ea) = (2.0, 35f64.to_radians(), 50.0, 3.0e4); + for kind in [ElementKind::Rod2, ElementKind::Rod3] { + let mut b = Build::statics(); + b.node(1, [-l * theta.tan(), l, 0.0]); + b.node(2, [0.0, l, 0.0]); + b.node(3, [l * theta.tan(), l, 0.0]); + b.node(4, [0.0, 0.0, 0.0]); + let mut members = Vec::new(); + for (i, top) in [1u32, 2, 3].into_iter().enumerate() { + let mut nodes = vec![top, 4]; + if kind == ElementKind::Rod3 { + // Middle nodes a third of the way along: still straight, not centred. + let (x0, x1) = (b.m.node(top).unwrap().xyz, [0.0; 3]); + let mid = 10 + i as u32; + b.node(mid, [0, 1, 2].map(|k| x0[k] + (x1[k] - x0[k]) / 3.0)); + nodes.push(mid); + // A straight member's middle node has no stiffness across it: hold it there. + let across = unit([x0[1], -x0[0], 0.0]); + let set = b.nodes(&format!("mid{mid}"), &[mid]); + b.spring(&format!("mid{mid}"), set, None, Direction::Vector(across), 1.0); + } + members.push((1 + i as u32, nodes)); + } + b.members("bars", kind, &members, 1.0, ea, None); + let top = b.nodes("top", &[1, 2, 3]); + b.fix("top", top); + let all: Vec = b.m.nodes().keys().copied().collect(); + let plane = b.nodes("plane_nodes", &all); + b.bc("plane", plane, ZERO_Z); + b.force("P", 4, [0.0, -p, 0.0]); + let r = b.solve(); + let c3 = theta.cos().powi(3); + let vertical = p / (1.0 + 2.0 * c3); + let inclined = p * theta.cos().powi(2) / (1.0 + 2.0 * c3); + let fe = member_forces(&b.m, &r, |_| 1.0); + let exact = [(1, inclined), (2, vertical), (3, inclined)]; + let err = worst_relative(&fe, &exact); + let delta = r.displacement[r.nodes.binary_search(&4).unwrap()][1]; + let want = -vertical * l / ea; + eprintln!("three-bar truss ({kind:?}): force error {err:.2e}, deflection {delta:.6e} vs {want:.6e}"); + assert!(err < 1e-10); + assert!((delta - want).abs() < 1e-10 * want.abs()); + } +} + +/// A rod3 bar with off-centre middle nodes under end tension carries σ = F/A exactly, at every +/// integration point and extrapolated to every node. +#[test] +fn rod3_constant_stress_patch() { + let (area, young, f) = (2.0, 1000.0, 6.0); + let mut b = Build::statics(); + let dir = unit([1.0, 2.0, 2.0]); + let at = |s: f64| dir.map(|c| c * s); + // Three elements of unequal length, middle nodes at 30%, 50% and 65% of each. + let ends = [0.0, 1.0, 2.5, 3.2]; + let fractions = [0.3, 0.5, 0.65]; + for (i, s) in ends.iter().enumerate() { + b.node(1 + i as u32, at(*s)); + } + let mut members = Vec::new(); + for i in 0..3 { + let s = ends[i] + fractions[i] * (ends[i + 1] - ends[i]); + b.node(10 + i as u32, at(s)); + members.push((1 + i as u32, vec![1 + i as u32, 2 + i as u32, 10 + i as u32])); + } + b.members("bar", ElementKind::Rod3, &members, area, young, None); + let start = b.nodes("start", &[1]); + b.fix("start", start); + // The other nodes slide along the bar only: hold them with springs across it. + let others: Vec = b.m.nodes().keys().copied().filter(|l| *l != 1).collect(); + let rest = b.nodes("rest", &others); + let across = unit([2.0, -1.0, 0.0]); + let third = [dir[1] * across[2] - dir[2] * across[1], dir[2] * across[0] - dir[0] * across[2], dir[0] * across[1] - dir[1] * across[0]]; + b.spring("across1", rest.clone(), None, Direction::Vector(across), 1.0e6); + b.spring("across2", rest, None, Direction::Vector(third), 1.0e6); + b.force("pull", 4, dir.map(|c| c * f)); + let r = b.solve(); + let mut worst: f64 = 0.0; + for e in &r.elements { + for p in &e.points { + worst = worst.max((axial(&p.stress, dir, false) - f / area).abs()); + } + for s in &e.nodal { + worst = worst.max((axial(s, dir, false) - f / area).abs()); + } + } + let tip = r.displacement[r.nodes.binary_search(&4).unwrap()]; + let stretch = tip.iter().zip(dir).map(|(u, d)| u * d).sum::(); + eprintln!("rod3 patch: stress error {worst:.2e}; tip stretch {stretch:.10} vs FL/EA {:.10}", f * 3.2 / (young * area)); + assert!(worst < 1e-10 * f / area); + assert!((stretch - f * 3.2 / (young * area)).abs() < 1e-12); +} + +/// A bar hanging under its own weight: u(z) = ρ g (L z − z²/2) / E from the support. Two-node +/// bars are exact at the nodes in one dimension; three-node bars reproduce the quadratic exactly. +#[test] +fn hanging_bar_under_gravity() { + let (l, rho, g, young, area) = (6.0, 2.0, 9.81, 1.0e4, 0.5); + for (kind, cells) in [(ElementKind::Rod2, 6usize), (ElementKind::Rod3, 3)] { + let mut b = Build::statics(); + let h = l / cells as f64; + for i in 0..=cells { + b.node(1 + i as u32, [0.0, 0.0, -(i as f64) * h]); + } + let mut members = Vec::new(); + for i in 0..cells { + let mut nodes = vec![1 + i as u32, 2 + i as u32]; + if kind == ElementKind::Rod3 { + b.node(100 + i as u32, [0.0, 0.0, -(i as f64 + 0.5) * h]); + nodes.push(100 + i as u32); + } + members.push((1 + i as u32, nodes)); + } + b.members("bar", kind, &members, area, young, Some(rho)); + let top = b.nodes("top", &[1]); + b.fix("top", top); + let all: Vec = b.m.nodes().keys().copied().collect(); + let line = b.nodes("line", &all); + b.bc("line", line, [Some(0.0), Some(0.0), None, None, None, None]); + b.load("weight", Load::Gravity { acceleration: [0.0, 0.0, -g], set: None }); + let r = b.solve(); + let mut worst: f64 = 0.0; + for (label, u) in r.nodes.iter().zip(&r.displacement) { + let z = -b.m.node(*label).unwrap().xyz[2]; + let exact = -rho * g * (l * z - z * z / 2.0) / young; + worst = worst.max((u[2] - exact).abs()); + } + let tip = rho * g * l * l / (2.0 * young); + eprintln!("hanging bar ({kind:?}): nodal displacement error {:.2e} of the tip's", worst / tip); + assert!(worst < 1e-12 * tip.max(1.0), "{worst}"); + assert!((r.reaction_total[2] - rho * g * area * l).abs() < 1e-10); + } +} + +/// One node, a spring to the ground in X and a mass: ω = √(k/m) (consistent and lumped mass +/// alike); a static force gives u = F/k. +#[test] +fn spring_mass_natural_frequency_and_static_stretch() { + let (k, m) = (400.0, 2.5); + for mass in [MassKind::Consistent, MassKind::Lumped] { + let mut b = Build::modes(1); + b.node(1, [1.0, 2.0, 3.0]); + let n = b.nodes("n", &[1]); + b.spring("k", n.clone(), None, Direction::Dof(1), k); + b.mass("m", n.clone(), Some(Mass::Isotropic(m)), None); + b.bc("yz", n, [None, Some(0.0), Some(0.0), None, None, None]); + let r = b.modes_of(mass); + let f = (k / m).sqrt() / (2.0 * PI); + let err = (r.modes[0].frequency - f).abs() / f; + eprintln!("spring–mass ({mass:?}): f {:.12} vs √(k/m)/2π {f:.12}, error {err:.1e}", r.modes[0].frequency); + assert!(err < 1e-12); + assert!((r.total_mass - m).abs() < 1e-12); + assert!((r.modes[0].effective_mass[0] - m).abs() < 1e-9); + } + let mut b = Build::statics(); + b.node(1, [0.0; 3]); + let n = b.nodes("n", &[1]); + b.spring("k", n.clone(), None, Direction::Dof(1), k); + b.bc("yz", n, [None, Some(0.0), Some(0.0), None, None, None]); + b.force("F", 1, [8.0, 0.0, 0.0]); + let r = b.solve(); + assert!((r.displacement[0][0] - 8.0 / k).abs() < 1e-15); + // The spring carries the load to the ground: no support reaction, and no imbalance. + assert_eq!(r.reaction_total[0], 0.0); + assert!(r.warnings.iter().all(|w| !w.contains("differ")), "{:?}", r.warnings); +} + +/// Springs to the ground along three orthogonal skew vectors: frequencies √(kᵢ/m) with the mode +/// shapes along the vectors. +#[test] +fn springs_in_skew_directions() { + let m = 3.0; + let a = unit([1.0, 1.0, 0.5]); + let b2 = unit([1.0, -1.0, 0.0]); + let c = [a[1] * b2[2] - a[2] * b2[1], a[2] * b2[0] - a[0] * b2[2], a[0] * b2[1] - a[1] * b2[0]]; + let ks = [30.0, 120.0, 270.0]; + let mut b = Build::modes(3); + b.node(7, [0.0; 3]); + let n = b.nodes("n", &[7]); + for (i, (v, k)) in [a, b2, c].iter().zip(ks).enumerate() { + b.spring(&format!("k{i}"), n.clone(), None, Direction::Vector(v.map(|x| 2.0 * x)), k); + } + b.mass("m", n, Some(Mass::Isotropic(m)), None); + let r = b.modes_of(MassKind::Consistent); + for ((mode, k), v) in r.modes.iter().zip(ks).zip([a, b2, c]) { + let f = (k / m).sqrt() / (2.0 * PI); + let s = mode.shape[0]; + let along = (s[0] * v[0] + s[1] * v[1] + s[2] * v[2]).abs() / (s[0] * s[0] + s[1] * s[1] + s[2] * s[2]).sqrt(); + eprintln!("skew springs: f {:.10} vs {f:.10}; shape along the spring {along:.12}", mode.frequency); + assert!((mode.frequency - f).abs() < 1e-10 * f); + assert!((along - 1.0).abs() < 1e-9); + } +} + +/// Ground — k₁ — m₁ — k₂ — m₂, the two-point spring acting along the line of its nodes: +/// ω² are the roots of m₁m₂ω⁴ − (k₁m₂ + k₂m₁ + k₂m₂)ω² + k₁k₂ = 0. +#[test] +fn two_mass_chain() { + let (k1, k2, m1, m2) = (100.0, 60.0, 2.0, 0.5); + let mut b = Build::modes(2); + let dir = unit([3.0, 0.0, 4.0]); + b.node(1, [0.0; 3]); + b.node(2, dir.map(|c| 2.0 * c)); + let n1 = b.nodes("n1", &[1]); + let n2 = b.nodes("n2", &[2]); + b.spring("k1", n1.clone(), None, Direction::Vector(dir), k1); + b.spring("k2", n1.clone(), Some(n2.clone()), Direction::Line, k2); + b.mass("m1", n1, Some(Mass::Isotropic(m1)), None); + b.mass("m2", n2, Some(Mass::Directional([m2, 7.0, m2])), None); + // Only the motion along the line: hold the others with stiff springs across it. + let both = b.nodes("both", &[1, 2]); + b.spring("hold_y", both.clone(), None, Direction::Dof(2), 1.0e9); + b.spring("hold_n", both, None, Direction::Vector([dir[2], 0.0, -dir[0]]), 1.0e9); + let r = b.modes_of(MassKind::Lumped); + let (pa, pb, pc) = (m1 * m2, -(k1 * m2 + k2 * m1 + k2 * m2), k1 * k2); + let disc = (pb * pb - 4.0 * pa * pc).sqrt(); + let exact = [(-pb - disc) / (2.0 * pa), (-pb + disc) / (2.0 * pa)]; + for (mode, w2) in r.modes.iter().zip(exact) { + let err = (mode.eigenvalue - w2).abs() / w2; + eprintln!("two-mass chain: ω² {:.10} vs {w2:.10}, error {err:.1e}", mode.eigenvalue); + assert!(err < 1e-8, "the holding springs perturb at 1e-9"); + } +} + +/// A massless bar with a tip mass M: ω = √(EA/(L M)); the bar's own mass ρAL adds the +/// Rayleigh-quotient correction ω² ≈ EA/(L(M + ρAL/3)) for a light bar. +#[test] +fn bar_with_a_tip_mass() { + let (l, young, area, tip) = (2.0, 2.0e5, 1e-2, 50.0); + for kind in [ElementKind::Rod2, ElementKind::Rod3] { + let mut b = Build::modes(1); + b.node(1, [0.0; 3]); + b.node(2, [l, 0.0, 0.0]); + let mut nodes = vec![1, 2]; + if kind == ElementKind::Rod3 { + b.node(3, [l / 2.0, 0.0, 0.0]); + nodes.push(3); + } + b.members("bar", kind, &[(1, nodes)], area, young, Some(1e-9)); + let root = b.nodes("root", &[1]); + b.fix("root", root); + let all: Vec = b.m.nodes().keys().copied().collect(); + let line = b.nodes("line", &all); + b.bc("line", line, [None, Some(0.0), Some(0.0), None, None, None]); + let end = b.nodes("end", &[2]); + b.mass("tip", end, Some(Mass::Isotropic(tip)), None); + let r = b.modes_of(MassKind::Consistent); + let f = (young * area / (l * tip)).sqrt() / (2.0 * PI); + let err = (r.modes[0].frequency - f).abs() / f; + eprintln!("bar with tip mass ({kind:?}): f {:.10} vs {f:.10}, error {err:.1e}", r.modes[0].frequency); + assert!(err < 1e-9); + } +} + +/// A torsional pendulum: a reference point with polar inertia J on a torsional spring k_t to +/// the ground; ω = √(k_t/J). The other directions are held. +#[test] +fn torsional_pendulum() { + let (kt, j) = (5.0e3, 0.8); + let mut b = Build::modes(1); + b.m.add_reference_point("Disc", ReferencePoint { xyz: [1.0, 1.0, 0.0], label: 0 }).unwrap(); + let disc = Target::ReferencePoint("Disc".into()); + b.spring("shaft", disc.clone(), None, Direction::Dof(6), kt); + // Products of inertia about the other axes do not couple into the held rotations. + b.mass("disc", disc.clone(), Some(Mass::Isotropic(4.0)), Some([0.5, 0.5, j, 0.0, 0.1, 0.1])); + b.bc("held", disc, [Some(0.0), Some(0.0), Some(0.0), Some(0.0), Some(0.0), None]); + let r = b.modes_of(MassKind::Lumped); + let f = (kt / j).sqrt() / (2.0 * PI); + let err = (r.modes[0].frequency - f).abs() / f; + eprintln!("torsional pendulum: f {:.12} vs √(k_t/J)/2π {f:.12}, error {err:.1e}", r.modes[0].frequency); + assert!(err < 1e-12); + assert_eq!(r.reference_points.len(), 1); + assert!(r.modes[0].point_shape[0][5].abs() > 0.0); +} + +/// Gravity acts on point masses (whole-model gravity); dashpots are ignored with a warning in +/// a static step; the runs are bit-for-bit repeatable. +#[test] +fn point_mass_gravity_dashpots_and_determinism() { + let mut b = Build::statics(); + b.node(1, [0.0; 3]); + b.node(2, [0.0, 0.0, -1.0]); + b.members("rod", ElementKind::Rod2, &[(1, vec![1, 2])], 1.0, 1.0e3, Some(0.0 + 1e-6)); + let top = b.nodes("top", &[1]); + b.fix("top", top); + let bottom = b.nodes("bottom", &[2]); + b.bc("line", bottom.clone(), [Some(0.0), Some(0.0), None, None, None, None]); + b.mass("weight", bottom.clone(), Some(Mass::Isotropic(3.0)), None); + b.m.add_spring("damper", Spring::Dashpot { first: bottom, second: None, direction: Direction::Dof(3), coefficient: 2.0 }).unwrap(); + b.load("g", Load::Gravity { acceleration: [0.0, 0.0, -10.0], set: None }); + let r = b.solve(); + let total = 3.0 * 10.0 + 1e-6 * 10.0; + assert!((r.reaction_total[2] - total).abs() < 1e-10, "{:?}", r.reaction_total); + assert!((r.displacement[r.nodes.binary_search(&2).unwrap()][2] + total * 1.0 / 1.0e3).abs() < 1e-6); + assert!(r.warnings.iter().any(|w| w.contains("dashpots `damper`")), "{:?}", r.warnings); + assert_eq!(r, b.solve()); +} + +/// What the solver cannot do is an error that says so. +#[test] +fn errors_name_the_problem() { + // A point mass with nothing holding it. + let mut b = Build::modes(1); + b.node(1, [0.0; 3]); + let n = b.nodes("n", &[1]); + b.mass("m", n, Some(Mass::Isotropic(1.0)), None); + // A free particle has no stiffness at all to shift against. + assert!(solve_frequency(&b.m, None, &FrequencyOptions::default(), &Limits::default()).is_err()); + // A truss of an orthotropic material. + let mut b = Build::statics(); + b.node(1, [0.0; 3]); + b.node(2, [1.0, 0.0, 0.0]); + b.m.set_element(1, ElementKind::Rod2, vec![1, 2]).unwrap(); + b.m.add_set("rod", Set { elements: vec![1], ..Set::default() }).unwrap(); + b.m.add_material( + "Ortho", + Material { + elastic: Some(Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [1.0; 3] }), + density: None, + ..Material::default() + }, + ) + .unwrap(); + b.m.add_section( + "rod", + Section { set: "rod".into(), material: "Ortho".into(), kind: SectionKind::Rod { area: 1.0 }, formulation: Formulation::default() }, + ) + .unwrap(); + let e = solve_static(&b.m, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("orthotropic"), "{e}"); +} + +/// A 2D (plane stress) plate stiffened by trusses along its loaded edges: the truss takes the +/// model's two in-plane unknowns per node, and the two act as parallel springs, +/// u = F L / (E t w + E_b A) (exact: uniform tension, ν = 0). +#[test] +fn plane_truss_stiffens_a_plane_stress_plate() { + let (l, w, t, young, area, young_bar, force) = (4.0, 1.0, 0.2, 1000.0, 0.05, 2.0e4, 30.0); + for (kind, mid) in [(ElementKind::Rod2, false), (ElementKind::Rod3, true)] { + let mut b = Build::statics(); + for (label, xyz) in [(1, [0.0, 0.0, 0.0]), (2, [l, 0.0, 0.0]), (3, [l, w, 0.0]), (4, [0.0, w, 0.0])] { + b.node(label, xyz); + } + b.m.set_element(1, ElementKind::Quad4, vec![1, 2, 3, 4]).unwrap(); + b.m.add_set("plate", Set { elements: vec![1], ..Set::default() }).unwrap(); + b.m.add_material("Plate", Material { elastic: Some(Elastic::Isotropic { young, poisson: 0.0 }), density: None, ..Material::default() }) + .unwrap(); + b.m.add_section( + "plate", + Section { + set: "plate".into(), + material: "Plate".into(), + kind: SectionKind::PlaneStress { thickness: t }, + formulation: Formulation::default(), + }, + ) + .unwrap(); + // Two bars of half the area, along the bottom and top edges, keep the tension uniform. + let mut members = vec![(2, vec![1, 2]), (3, vec![4, 3])]; + if mid { + b.node(5, [l / 2.0, 0.0, 0.0]); + b.node(6, [l / 2.0, w, 0.0]); + members[0].1.push(5); + members[1].1.push(6); + } + b.members("bar", kind, &members, area / 2.0, young_bar, None); + let left = b.nodes("left", &[1, 4]); + b.bc("left", left, [Some(0.0), None, None, None, None, None]); + let corner = b.nodes("corner", &[1]); + b.bc("corner", corner, [None, Some(0.0), None, None, None, None]); + if mid { + let middle = b.nodes("middle", &[5, 6]); + b.bc("middle", middle, [None, Some(0.0), None, None, None, None]); + } + let right = b.nodes("right", &[2, 3]); + b.load("pull", Load::Force { target: right, vector: [force, 0.0, 0.0], distribution: Distribution::PerNode }); + let r = b.solve(); + let u = r.displacement[r.nodes.binary_search(&2).unwrap()][0]; + let exact = 2.0 * force * l / (young * t * w + young_bar * area); + eprintln!("plane truss ({kind:?}): u {u:.12} vs {exact:.12}; {} unknowns", r.stats.dofs); + assert!((u - exact).abs() < 1e-12 * exact); + assert_eq!(r.stats.dofs, 2 * b.m.nodes().len(), "two unknowns per node in a 2D model"); + } +} diff --git a/crates/ui-egui/src/theme.rs b/crates/ui-egui/src/theme.rs index d56f8f0..55dda3e 100644 --- a/crates/ui-egui/src/theme.rs +++ b/crates/ui-egui/src/theme.rs @@ -153,7 +153,7 @@ impl Tokens { ElementKind::Tet4 | ElementKind::Tet10 => self.tet, ElementKind::Wedge6 => self.wedge, ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => self.shell, - ElementKind::Rod2 | ElementKind::Beam2 => self.line_element, + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => self.line_element, } } } diff --git a/docs/architecture.md b/docs/architecture.md index 80c0278..ece476f 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -138,8 +138,13 @@ module docs; ADR 0001 records why the linear algebra is our own. quad4, quad8), chosen by the section kind, and the model then has two degrees of freedom per node: the 2D elements declare `ElementDofs::new(DofSet::PLANE)` (Ux, Uy) to the `DofLayout`. A 2D model lies in the XY plane (axisymmetric: x the radius, Y the axis; nodes on the axis are held - radially), is loaded on element edges, and cannot mix with solids or couplings yet. Every reference point a coupling - ties to the mesh gets six degrees of freedom (translations, rotations) after the mesh's. Shells + radially), is loaded on element edges, and cannot mix with solids or couplings yet. Rod + sections' rod2 and rod3 elements are trusses (`modelcraft_fem::truss`). Point masses, rotary + inertia and springs (`Resolved::masses`, `springs`, expanded per node by the model) assemble on + the unknowns their nodes are given (`discrete::needs`, `discrete::dof_of`); dashpots are + ignored with a warning (no velocities in static and frequency steps). Every reference point a + coupling or a discrete element uses gets six degrees of freedom (translations, rotations) after + the mesh's. Shells and beams return "not supported yet". The step's boundary conditions prescribe translations (rotations do not exist on solid nodes) and a coupled point's six directions; its loads become a load vector: forces and tractions from `Resolved::nodal_loads`, forces and moments on coupled @@ -181,7 +186,10 @@ and rigid-body modes. (kinematic and distributing), clamping through a coupled point equals clamping (static and modes), tension through a distributing coupling against the equivalent pressure away from the end, pure bending by a moment through a distributing coupling, torsion of a tube by a moment on -a kinematically coupled point (twist and shear stress). +a kinematically coupled point (twist and shear stress). `crates/solver/tests/one_d.rs`: member +forces of determinate trusses (a plane Pratt truss, a space tripod) against the method of joints, +the three-bar truss, a rod3 constant-stress patch and a bar hanging under its weight, spring–mass, +skew springs, the two-mass chain, a bar with a tip mass and the torsional pendulum. ## Explicit solver @@ -192,9 +200,11 @@ number their degrees of freedom with `modelcraft_fem::DofLayout` (docs/fem.md). 1. **System.** `Model::resolve()`; hex8 sections with reduced integration become one-point uniform-strain hexahedra with Flanagan–Belytschko hourglass control, tet4 sections Bonet–Burton - average nodal pressure tetrahedra (`modelcraft_fem::uniform` holds their kinematics). Other - element kinds and formulations, couplings, reference points and moments return "not supported - yet". Masses are row-sum lumped. A zero-valued boundary condition holds a degree of freedom; a + average nodal pressure tetrahedra (`modelcraft_fem::uniform` holds their kinematics), rod2 and + rod3 trusses (`discrete`). Point masses add to the lumped mass; springs and dashpots (on the + previous half increment's velocity) add internal and viscous forces and their own stable + increment. Other element kinds and formulations, couplings, reference points, rotary inertia, + rotational springs and moments return "not supported yet". Masses are row-sum lumped. A zero-valued boundary condition holds a degree of freedom; a displacement value is reached at constant velocity over the step. Loads act in full from the start. The step lists initial velocities and fixed frictionless rigid planes by set. 2. **Time step.** Each element's highest frequency is bounded in closed form (exact for a cube); @@ -227,7 +237,9 @@ only): stress waves in hex8 and tet4 bars (front at c·t, σ = ρ c v, work inpu rigid wall (contact time 2L/c, rebound speed, plateau force ρ c v₀ A), the energy balance, bit-for-bit repeatability; Dienes' simple shear for both rates, a stressed element turned rigidly by 90°, the J2 uniaxial curve and unloading at large strain, damage failure at the fracture -strain, and the Taylor bar against Taylor's 1948 estimate; in the engine (`solve.rs`), a +strain, and the Taylor bar against Taylor's 1948 estimate; `crates/explicit/tests/one_d.rs`: +spring–mass period, damped decay rate and period, a dashpot alone, two masses on a line spring, +waves in rod2 and rod3 bars; in the engine (`solve.rs`), a cantilever's free-vibration period against `solve.frequency`. ## Lattice Boltzmann solver diff --git a/docs/elements.md b/docs/elements.md index 0ede686..9a59649 100644 --- a/docs/elements.md +++ b/docs/elements.md @@ -19,10 +19,10 @@ hex8/20/27 with first and second derivatives). | Element | Use | Status | |---|---|---| -| Point mass | Lumped mass at a node or reference point (isotropic or per direction) | in progress (1D elements I) | -| Rotary inertia | Inertia tensor at a 6-DOF node | in progress (1D elements I) | -| Spring to ground / two-point spring | Linear spring in a fixed direction or along the line of nodes | in progress (1D elements I) | -| Dashpot to ground / two-point | Linear viscous damper (dynamic steps; ignored in static) | in progress (1D elements I) | +| Point mass | Lumped mass at a node or reference point (isotropic or per direction) | done (implicit, explicit; Abaqus MASS, Nastran CONM1, LS-DYNA ELEMENT_MASS; per direction: Nastran only) | +| Rotary inertia | Inertia tensor at a 6-DOF node | done (implicit; Abaqus ROTARYI, Nastran CONM1, LS-DYNA ELEMENT_INERTIA; explicit when it has rotations) | +| Spring to ground / two-point spring | Linear spring in a fixed direction or along the line of nodes | done (implicit, explicit translational; Abaqus SPRING1/2/A, Nastran CELAS2, LS-DYNA ELEMENT_DISCRETE; skew fixed directions: LS-DYNA only) | +| Dashpot to ground / two-point | Linear viscous damper (dynamic steps; ignored in static) | done (explicit; ignored with a warning in static and frequency steps; Abaqus DASHPOT1/2/A, Nastran CDAMP2, LS-DYNA ELEMENT_DISCRETE) | | Nonstructural mass | Mass spread over a region without stiffness | planned | | Connectors | Kinematic joints (hinge, slot, universal…) with elasticity, damping, stops | planned (M7) | @@ -30,8 +30,8 @@ hex8/20/27 with first and second derivatives). | Element | Use | Status | |---|---|---| -| Truss, 2-node (2D and 3D) | Axial members, cables (compression allowed) | in progress (1D elements I) | -| Truss, 3-node (2D and 3D) | Curved axial members | in progress (1D elements I) | +| Truss, 2-node (2D and 3D) | Axial members, cables (compression allowed) | done (rod2: implicit, in plane models with their two unknowns per node; explicit, corotational under nonlinear geometry; Abaqus T3D2 / T2D2, Nastran CROD, LS-DYNA truss beam (3D models)) | +| Truss, 3-node (2D and 3D) | Curved axial members | done (rod3: implicit, plane models too; explicit with small displacements; Abaqus T3D3 / T2D3; no Nastran or LS-DYNA counterpart) | | Beam, 2-node Timoshenko, 3D | Frames, stiffeners; shear deformable | in progress (shells and beams) | | Beam, 2-node Euler–Bernoulli (cubic), 2D and 3D | Slender frames; exact for end-loaded prismatic members | planned (1D elements II) | | Beam, 3-node Timoshenko (quadratic), 2D and 3D | Curved and stocky beams | planned (1D elements II) | diff --git a/docs/features.md b/docs/features.md index 85421da..ac5ce75 100644 --- a/docs/features.md +++ b/docs/features.md @@ -5,7 +5,7 @@ Generated from the capability catalogue (`crates/engine/src/catalogue.rs`) by Every entry is a toolbar button. **Ready** entries run their command; **coming soon** entries are shown greyed with the roadmap milestone that brings them (ROADMAP.md). -48 ready, 53 coming soon. +50 ready, 51 coming soon. ## Geometry @@ -106,8 +106,8 @@ shown greyed with the roadmap milestone that brings them (ROADMAP.md). | Contact | Find Contacts — Detect contact pairs automatically | coming soon (M7) | | | Contact | Rigid Wall — Rigid walls for explicit dynamics | coming soon (M10) | | | Connectors | Bolt Load — Bolt pretension | coming soon (M7) | | -| Connectors | Spring — Springs and dashpots | coming soon (M7) | | -| Connectors | Point Mass — Point masses and inertia | coming soon (M7) | | +| Connectors | Spring — Spring or dashpot to the ground or between two points, along a direction or their line | ready | `model.spring.create`, `model.spring.edit`, `model.spring.delete` | +| Connectors | Point Mass — Point mass and rotary inertia at nodes or a reference point | ready | `model.mass.create`, `model.mass.edit`, `model.mass.delete` | | Connectors | Connector — Connector elements between points | coming soon (M7) | | | Coupled | FSI Interface — The fluid–structure interface | coming soon (M13) | | diff --git a/docs/fem.md b/docs/fem.md index 38d0a86..3e2fb41 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -39,6 +39,12 @@ Lookups: `node(label)` / `slot(label)` (a slot is a position: mesh nodes `0..nod `point_slot(k)`), `dof(slot, kind)`, `range(slot)`, `kinds(slot)`, `locate(dof) -> (slot, kind)`, `element_dofs(slots, kinds)`, `describe(dof)`, `stride()` (uniform layouts). +Trusses declare `modelcraft_fem::truss::DOFS` (translations, on `Basis::Line2`/`Line3`). +Point masses, springs and dashpots are not mesh elements: the implicit solver asks the layout +for what each needs at its nodes (`modelcraft_solver::discrete::needs`: translations for a mass +or a translational spring, rotations for an inertia or a torsional spring; a reference point +they use is a point with six) and finds its unknowns with `discrete::dof_of`. + Solvers: `modelcraft_solver::system::Mesh::layout` numbers the implicit system (mesh nodes, then six per coupled reference point); assembly (`sparse::add_block(k, layout, slots, kinds, data)`, `add_scalar`, `add_diagonal`) works for any mix of kinds per node, so a solid sharing diff --git a/docs/parity.md b/docs/parity.md index 4352e52..54cd3a6 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -56,7 +56,8 @@ products are run, no documentation text is copied). Status: **done**, **started* | Model | Steps: nonlinear static, explicit dynamics, heat transfer; amplitudes | M7–M8 | | planned | | Interaction | Tie between non-matching meshes | M7 | | planned | | Interaction | Contact (general and pairs), friction, automatic pair detection | M7 | | planned | -| Interaction | Bolt pretension, springs, point masses, connectors | M7 | | planned | +| Interaction | Springs and dashpots (ground or two-point; DOF, vector or line), point masses with rotary inertia | M4 | | done (implicit; explicit translational; all three decks with what each cannot take reported) | +| Interaction | Bolt pretension, connectors | M7 | | planned | | Decks | Abaqus input writer (`deck.write`, flat deck; reports what Abaqus cannot express) | M3 | | done | | Decks | Nastran bulk data writer (`deck.write`, SOL 101 or 103 per deck; reports what Nastran cannot express) | M3 | | done | | Decks | LS-DYNA keyword writer (implicit; explicit in M8) | M3 | | done (implicit static and eigenvalue; beam sections planned) | diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index 68a3537..f87fec7 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -997,6 +997,12 @@ the module docs): distributing (interpolation) constraint as a least-squares rigid fit with loads spread by virtual work, G. H. K. Heirman, W. Desmet, Multibody Syst. Dyn. 24(2):219–234, 2010; the reduced matrix product T_fᵀ K T_f row by row, F. G. Gustavson, ACM TOMS 4(3):250–269, 1978. +- Trusses and discrete elements: the isoparametric bar along its arc, R. D. Cook, D. S. Malkus, + M. E. Plesha, R. J. Witt, *Concepts and Applications of Finite Element Analysis*, 4th ed., + Wiley 2002, §3.6 and §6.2, and Hughes 2000, §1.10; springs k [n nᵀ, −n nᵀ; −n nᵀ, n nᵀ], + Cook et al. 2002, §2.2 and §2.10; HRZ lumping as above. Verified against the method of joints + (statically determinate trusses), the three-bar truss, the hanging bar, √(k/m), the two-mass + chain and the torsional pendulum √(k_t/J). - Verification: R. H. MacNeal, R. L. Harder, *A proposed standard set of problems to test finite element accuracy*, Finite Elem. Anal. Des. 1(1):3–20, 1985 (patch tests); Lamé, Kirsch, Timoshenko and Euler–Bernoulli closed forms; Saint-Venant's principle (tension through a @@ -1024,6 +1030,10 @@ again in the module docs): 14(5):437–449, 1998. - Bulk viscosity and its time-step correction: J. von Neumann, R. D. Richtmyer, J. Appl. Phys. 21(3):232–237, 1950; M. L. Wilkins, J. Comput. Phys. 36(3):281–303, 1980. +- Springs and dashpots: lagged-velocity damping and the damped stable increment + Δt ≤ (2/ω)(√(1 + ξ²) − ξ), Belytschko et al. 2014, §6.2–6.3; Gershgorin's bound on + M^{-1/2} K M^{-1/2} per degree of freedom. Verified against the damped oscillator's decay e^{−ξωt} + and period 2π/ω_d. - Rigid walls as kinematic velocity constraints: T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier, W. Kanoknukulchai, CMAME 8(3):249–276, 1976. - Verification: d'Alembert's solution of the 1D wave equation (front at c t, σ = ρ c v), the