diff --git a/Cargo.lock b/Cargo.lock index 748d643..17865ae 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1846,6 +1846,7 @@ name = "modelcraft-decks" version = "0.1.0" dependencies = [ "modelcraft-model", + "serde_json", "thiserror 2.0.21", ] @@ -1854,6 +1855,7 @@ name = "modelcraft-engine" version = "0.1.0" dependencies = [ "log", + "modelcraft-decks", "modelcraft-geometry", "modelcraft-io", "modelcraft-mesh", diff --git a/ROADMAP.md b/ROADMAP.md index 8e1184e..c8b0e21 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -21,7 +21,7 @@ hex meshes today and it gives the best elements where stresses are read. | M0 Foundation | Workspace, quality gates, FE model, command engine with rollback and panic guard, CLI, geometry boundary on SolidCraft | done | | 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); 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, Nastran SOL 101/103 and LS-DYNA implicit; 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 | +| 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 and LS-DYNA implicit; 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 | Linear static and modal solver (tet4/10, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions; deformed shape, contours, mode shapes, probes, reaction sums | planned | | 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 | diff --git a/crates/decks/Cargo.toml b/crates/decks/Cargo.toml index 97bf2c1..cc0d16e 100644 --- a/crates/decks/Cargo.toml +++ b/crates/decks/Cargo.toml @@ -11,3 +11,6 @@ workspace = true [dependencies] modelcraft-model = { workspace = true } thiserror = { workspace = true } + +[dev-dependencies] +serde_json = { workspace = true } diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs new file mode 100644 index 0000000..f38fc8a --- /dev/null +++ b/crates/decks/src/abaqus.rs @@ -0,0 +1,796 @@ +//! Abaqus input (`.inp`) writer. +//! +//! The deck is written from the public Abaqus keyword and element references, in our own words; +//! nothing here comes from a deck produced by the product. +//! +//! # Layout +//! +//! The deck is *flat*: nodes, elements, sets, materials, sections and couplings sit at model level, +//! with no parts, assembly or instances. ModelCraft has one mesh with model-wide labels and no +//! instance transforms, so a part/assembly wrapper would only prefix every set name with an +//! instance name (`Part-1-1.root`) and bring the 80-character limit on the joined name closer. When +//! the model gains several parts or instances the writer switches to the part/assembly form. +//! +//! In order: heading, nodes (mesh nodes, then reference points), elements grouped by element type, +//! sets (a node set, an element set and an element-face surface per model set, each written when +//! the set resolves to some), materials, sections, couplings, the initial state's boundary +//! conditions, then the steps. +//! +//! # Steps and the activity table +//! +//! The model lists, for every step, exactly the boundary conditions and loads active in it; Abaqus +//! carries them from step to step. The writer bridges the two: +//! - a boundary condition active in every step that holds its degrees of freedom at zero is written +//! once, in the initial state (model data); +//! - a general (static) step whose active boundary conditions equal the previous general step's +//! writes none; one that only adds writes the additions (`OP=MOD`); otherwise it writes +//! `*BOUNDARY, OP=NEW` and restates every active one; +//! - loads in a general step are restated, keyword by keyword, with `OP=NEW` whenever the active set +//! changes, so a load not listed is removed; +//! - a frequency step is a perturbation of the last general step: its boundary conditions are +//! written against that base state and do not carry on to later steps; loads have no effect on +//! natural frequencies and are not written in it; +//! - every step restates its output request (`*OUTPUT, FIELD, OP=NEW`). +//! +//! # Loads +//! +//! - Forces, moments and tractions are written as concentrated loads from the model's resolved +//! nodal values (the same distribution every writer and the solver use), summed per node and +//! degree of freedom so loads sharing a node add up. +//! - Pressure is a surface load on the set's element-face surface, summed per set. +//! - Gravity and body forces are element loads on element sets (the section sets when gravity acts +//! on the whole model); when two active loads of one kind share elements they are summed and +//! written per element. +//! +//! # Numbers and labels +//! +//! Real numbers are written in the shortest form that reads back to the same `f64` (exponent marker +//! `E`). Lines stay within 256 characters: element connectivity puts at most 15 nodes after the +//! label on the first line and 16 per continuation line; sets list 16 labels per line. Labels are +//! the model's names (a letter, then letters, digits, `_` or `-`, at most 64 characters), which +//! Abaqus reads unquoted. Abaqus ignores case in names, so two names that differ only in case are +//! reported as [`Error::Unsupported`], as is everything else the format cannot express. +//! +//! # Element types +//! +//! | ModelCraft | Solid section | Shell section | Beam | Rod | +//! |---|---|---|---|---| +//! | tet4 | C3D4 (H) | | | | +//! | tet10 | C3D10 (H) | | | | +//! | wedge6 | C3D6 (H) | | | | +//! | hex8 | C3D8, C3D8R, C3D8I (H) | | | | +//! | hex20 | C3D20, C3D20R (H) | | | | +//! | tri3, quad4 | | S3 / S3R, S4 / S4R | | | +//! | quad8 | | S8R (reduced only) | | | +//! | beam2, rod2 | | | B31 | T3D2 | +//! +//! `H` marks the hybrid variant. A four-node tetrahedron integrates at one point either way, so +//! reduced integration writes C3D4. Not expressible: reduced integration of ten-node tetrahedra and +//! six-node wedges, full integration of eight-node shells, six-node triangular shells, orthotropic +//! materials outside solid sections (shell and beam material axes need orientations), and +//! distributing couplings that leave a translation free. +//! +//! Node order: the model stores elements in Abaqus order, so connectivity passes through; the +//! model's face numbers are Abaqus's S1–S6 for solids, and SPOS/SNEG for shell faces 1/2. + +use std::collections::{BTreeMap, BTreeSet}; +use std::fmt::Write as _; + +use modelcraft_model::{ + BeamProfile, BoundaryCondition, CouplingKind, Elastic, Element, ElementKind, Integration, Load, Model, OutputVariable, Procedure, Resolved, + Section, SectionKind, Step, Symmetry, Target, +}; + +use crate::{Error, Result, unsupported}; + +const FORMAT: &str = "Abaqus"; + +/// Longest line Abaqus reads. +pub const MAX_LINE: usize = 256; + +/// Labels per set data line. +const SET_LABELS_PER_LINE: usize = 16; + +/// Nodes on an element's first data line (after its label), and on each continuation line. +const FIRST_LINE_NODES: usize = 15; +const CONTINUATION_NODES: usize = 16; + +/// Write the model as an Abaqus input deck. Fails when the model does not resolve (its check +/// finds an error) or has something the format cannot express ([`Error::Unsupported`], every such +/// item named). +pub fn write(model: &Model) -> Result { + write_resolved(model, &model.resolve()?) +} + +/// [`write`] for a model already resolved. +pub fn write_resolved(model: &Model, r: &Resolved) -> Result { + let by_type = elements_by_type(model, r)?; + check_expressible(model)?; + let mut d = Deck::default(); + d.heading(model); + d.nodes(model)?; + d.elements(&by_type); + d.sets(model, r)?; + d.materials(model)?; + d.sections(model)?; + d.couplings(model, r); + d.steps(model, r)?; + Ok(d.out) +} + +/// The Abaqus element type of an element of `kind` in `section`, or why there is none. +pub fn element_type(kind: ElementKind, section: &Section) -> std::result::Result { + use ElementKind as K; + use SectionKind as S; + let f = section.formulation; + let reduced = f.integration == Integration::Reduced; + if (f.hybrid || f.incompatible_modes) && section.kind.dimension() != 3 { + return Err(format!("hybrid or incompatible-mode {} sections", section.kind.name())); + } + let plain = match (§ion.kind, kind, reduced, f.incompatible_modes) { + (S::Solid {}, K::Hex8, false, true) => "C3D8I", + (S::Solid {}, _, true, true) => return Err("incompatible modes with reduced integration".into()), + (S::Solid {}, _, _, true) => return Err(format!("incompatible modes on {} elements (hex8 only)", kind.name())), + (S::Solid {}, K::Tet4, ..) => "C3D4", + (S::Solid {}, K::Tet10 | K::Wedge6, true, _) => { + return Err(format!("reduced integration of {} elements (Abaqus has no reduced-integration variant)", kind.name())); + } + (S::Solid {}, K::Tet10, ..) => "C3D10", + (S::Solid {}, K::Wedge6, ..) => "C3D6", + (S::Solid {}, K::Hex8, r, _) => ["C3D8", "C3D8R"][usize::from(r)], + (S::Solid {}, K::Hex20, r, _) => ["C3D20", "C3D20R"][usize::from(r)], + (S::Shell { .. }, K::Tri3, r, _) => ["S3", "S3R"][usize::from(r)], + (S::Shell { .. }, K::Quad4, r, _) => ["S4", "S4R"][usize::from(r)], + (S::Shell { .. }, K::Quad8, true, _) => "S8R", + (S::Shell { .. }, K::Quad8, false, _) => { + return Err("full integration of quad8 shells (Abaqus's eight-node shell, S8R, uses reduced integration)".into()); + } + (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, ..) => return Err(format!("{} elements in a {} section", kind.name(), s.name())), + }; + Ok(if f.hybrid { format!("{plain}H") } else { plain.to_string() }) +} + +/// Every element, grouped by element type (from its section) and sorted by label. +fn elements_by_type<'m>(model: &'m Model, r: &Resolved) -> Result>> { + let mut by_type: BTreeMap> = BTreeMap::new(); + let mut bad = BTreeSet::new(); + let mut count = 0; + for (name, section) in model.sections() { + let mut by_kind: BTreeMap> = BTreeMap::new(); + for &e in r.sets.get(§ion.set).map(|s| s.elements.as_slice()).unwrap_or_default() { + if let Some(element) = model.element(e) { + by_kind.entry(element.kind).or_default().push((e, element)); + } + } + for (kind, elements) in by_kind { + match element_type(kind, section) { + Ok(t) => { + count += elements.len(); + by_type.entry(t).or_default().extend(elements); + } + Err(why) => { + bad.insert(format!("section `{name}`: {why}")); + } + } + } + } + unsupported(FORMAT, bad)?; + // A resolved model puts every element in exactly one section; this guards the writer. + if count != model.elements().len() { + return Err(Error::Unsupported { format: FORMAT, what: "elements without exactly one section".into() }); + } + for elements in by_type.values_mut() { + elements.sort_unstable_by_key(|(e, _)| *e); + } + Ok(by_type) +} + +/// Everything else the format cannot express: case-only name clashes, orthotropic materials +/// outside solids, distributing couplings with a free translation. +fn check_expressible(model: &Model) -> Result<()> { + let mut bad = BTreeSet::new(); + let namespaces: [(&str, Vec<&str>); 4] = [ + ("sets", model.sets().keys().map(String::as_str).collect()), + ("materials", model.materials().keys().map(String::as_str).collect()), + ("couplings", model.couplings().keys().map(String::as_str).collect()), + ("steps", model.steps().iter().map(|s| s.name.as_str()).collect()), + ]; + for (what, names) in namespaces { + let mut seen: BTreeMap = BTreeMap::new(); + for n in names { + if let Some(other) = seen.insert(n.to_ascii_uppercase(), n) { + bad.insert(format!("{what} `{other}` and `{n}` differ only in case, and Abaqus ignores case in names")); + } + } + } + for (name, section) in model.sections() { + let ortho = matches!(model.material(§ion.material).and_then(|m| m.elastic), Some(Elastic::EngineeringConstants { .. })); + if ortho && section.kind.dimension() != 3 { + bad.insert(format!( + "section `{name}`: orthotropic material `{}` in a {} section (its material axes need an orientation, not supported yet)", + section.material, + section.kind.name() + )); + } + } + for (name, c) in model.couplings() { + if c.kind == CouplingKind::Distributing && !(1..=3).all(|d| c.dofs.contains(d)) { + bad.insert(format!("coupling `{name}`: an Abaqus distributing coupling always couples translations 1-3")); + } + } + unsupported(FORMAT, bad) +} + +/// 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<()> { + if !v.is_finite() { + return Err(Error::Unsupported { format: FORMAT, what: format!("the value {v} (not a finite number)") }); + } + // Rust's `Debug` for f64 is the shortest round-trip form, switching to an exponent for very + // large and small magnitudes; Abaqus wants the marker `E`, and a mantissa with a point reads + // the same everywhere. + let start = out.len(); + let _ = write!(out, "{v:?}"); + if let Some(at) = out.get(start..).and_then(|s| s.find('e')).map(|i| start + i) { + let point = out.get(start..at).is_some_and(|m| m.contains('.')); + out.replace_range(at..at + 1, if point { "E" } else { ".0E" }); + } + Ok(()) +} + +/// [`push_number`] into a new string. +pub fn number(v: f64) -> Result { + let mut s = String::new(); + push_number(&mut s, v)?; + Ok(s) +} + +/// The surface face identifier of face `face` (1-based) of an element of `kind`. +pub fn face_id(kind: ElementKind, face: u8) -> Option<&'static str> { + const SOLID: [&str; 6] = ["S1", "S2", "S3", "S4", "S5", "S6"]; + const SHELL: [&str; 2] = ["SPOS", "SNEG"]; + let i = usize::from(face).checked_sub(1)?; + if i >= kind.faces().len() { + return None; + } + match kind.dimension() { + 3 => SOLID.get(i).copied(), + 2 => SHELL.get(i).copied(), + _ => None, + } +} + +/// The node surface each coupling ties its reference point to: a fresh name per coupling that +/// clashes with no set (sets with faces own surfaces of their name) and no other coupling surface. +fn coupling_surfaces(model: &Model, r: &Resolved) -> BTreeMap { + let mut taken: BTreeSet = r.sets.iter().filter(|(_, s)| !s.faces.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 +} + +/// Consecutive runs of degrees of freedom: `[1, 2, 3, 5]` → `[(1, 3), (5, 5)]`. +fn dof_runs(dofs: &[u8]) -> Vec<(u8, u8)> { + let mut runs: Vec<(u8, u8)> = Vec::new(); + for &d in dofs { + match runs.last_mut() { + Some((_, last)) if *last + 1 == d => *last = d, + _ => runs.push((d, d)), + } + } + runs +} + +/// Whether a boundary condition holds every degree of freedom it constrains at zero. +fn holds_at_zero(bc: &BoundaryCondition) -> bool { + (1..=6).all(|d| bc.value(d).is_none_or(|v| v == 0.0)) +} + +/// Where an element load acts: an element set, or one element. +#[derive(Clone, Copy)] +enum ElementTarget<'a> { + Set(&'a str), + Element(u32), +} + +/// Gravity (an acceleration) or a body force (a force per volume). +#[derive(Clone, Copy, PartialEq)] +enum BodyKind { + Gravity, + Force, +} + +/// A gravity or body force on some elements, before it is written. +struct ElementLoad<'a> { + /// The element sets that hold exactly its elements. + elsets: Vec<&'a str>, + elements: &'a [u32], + vector: [f64; 3], +} + +impl<'a> ElementLoad<'a> { + fn on_set(r: &'a Resolved, set: &'a str, vector: [f64; 3]) -> ElementLoad<'a> { + ElementLoad { elsets: vec![set], elements: r.sets.get(set).map(|s| s.elements.as_slice()).unwrap_or_default(), vector } + } +} + +/// Whether no element is in two of the loads. +fn disjoint(loads: &[ElementLoad]) -> bool { + if loads.len() < 2 { + return true; + } + let mut all: Vec = loads.iter().flat_map(|l| l.elements.iter().copied()).collect(); + let n = all.len(); + all.sort_unstable(); + all.dedup(); + all.len() == n +} + +/// Element loads of one kind: by element set when the loads share no element, otherwise summed +/// per element. Zero vectors are left out. +fn element_loads<'a>(loads: &[ElementLoad<'a>], kind: BodyKind, out: &mut Vec<(ElementTarget<'a>, BodyKind, [f64; 3])>) { + if disjoint(loads) { + for l in loads { + out.extend(l.elsets.iter().map(|s| (ElementTarget::Set(s), kind, l.vector))); + } + } else { + let mut summed: BTreeMap = BTreeMap::new(); + for l in loads { + for &e in l.elements { + for (c, add) in summed.entry(e).or_default().iter_mut().zip(l.vector) { + *c += add; + } + } + } + out.extend(summed.into_iter().map(|(e, v)| (ElementTarget::Element(e), kind, v))); + } + out.retain(|(_, _, v)| v.iter().any(|c| *c != 0.0)); +} + +/// The loads of one general step, as the deck writes them. +#[derive(Default)] +struct StepLoads<'a> { + /// (node, degree of freedom) → summed value. + concentrated: BTreeMap<(u32, u8), f64>, + /// Surface (set name) → summed pressure. + pressure: BTreeMap<&'a str, f64>, + element: Vec<(ElementTarget<'a>, BodyKind, [f64; 3])>, +} + +/// Which load keywords a step wrote entries for, so the next step knows what to remove. +#[derive(Default, Clone, Copy)] +struct Written { + concentrated: bool, + pressure: bool, + element: bool, +} + +impl<'a> StepLoads<'a> { + /// `everything` is every element and the section sets that hold them (gravity on the whole + /// model). + fn new(model: &'a Model, r: &'a Resolved, step: &'a Step, everything: &'a (Vec<&'a str>, Vec)) -> StepLoads<'a> { + let mut s = StepLoads { concentrated: crate::concentrated_loads(model, r, &step.loads), ..StepLoads::default() }; + let (mut gravity, mut body) = (Vec::new(), Vec::new()); + for name in &step.loads { + match model.load(name) { + Some(Load::Pressure { set, magnitude }) => *s.pressure.entry(set.as_str()).or_default() += magnitude, + Some(Load::Gravity { acceleration, set: Some(set) }) => gravity.push(ElementLoad::on_set(r, set, *acceleration)), + Some(Load::Gravity { acceleration, set: None }) => { + gravity.push(ElementLoad { elsets: everything.0.clone(), elements: &everything.1, vector: *acceleration }) + } + Some(Load::BodyForce { set, vector }) => body.push(ElementLoad::on_set(r, set, *vector)), + _ => {} + } + } + s.pressure.retain(|_, p| *p != 0.0); + element_loads(&gravity, BodyKind::Gravity, &mut s.element); + element_loads(&body, BodyKind::Force, &mut s.element); + s + } + + fn written(&self) -> Written { + Written { concentrated: !self.concentrated.is_empty(), pressure: !self.pressure.is_empty(), element: !self.element.is_empty() } + } +} + +#[derive(Default)] +struct Deck { + out: String, +} + +impl Deck { + fn line(&mut self, s: &str) { + self.out.push_str(s); + self.out.push('\n'); + } + + fn comment(&mut self, s: &str) { + self.out.push_str("** "); + self.line(s); + } + + /// `, value` for each value, then the end of the line. + fn numbers_line(&mut self, head: &str, values: &[f64]) -> Result<()> { + self.out.push_str(head); + for (i, v) in values.iter().enumerate() { + if i > 0 || !head.is_empty() { + self.out.push_str(", "); + } + push_number(&mut self.out, *v)?; + } + self.out.push('\n'); + Ok(()) + } + + /// Labels, [`SET_LABELS_PER_LINE`] per line. + fn labels(&mut self, labels: &[u32]) { + for chunk in labels.chunks(SET_LABELS_PER_LINE) { + for (i, l) in chunk.iter().enumerate() { + let _ = write!(self.out, "{}{l}", if i > 0 { ", " } else { "" }); + } + self.out.push('\n'); + } + } + + fn heading(&mut self, model: &Model) { + let units = match model.units().length() { + Some(_) => format!("units {}", model.units().label()), + None => "units not specified".to_string(), + }; + self.line("*HEADING"); + self.line(&format!("ModelCraft model, {units}")); + self.comment("Written by ModelCraft. The numbers are in one consistent unit system and are not converted."); + } + + fn nodes(&mut self, model: &Model) -> Result<()> { + self.comment("Nodes"); + self.line("*NODE"); + for (label, node) in model.nodes() { + self.numbers_line(&label.to_string(), &node.xyz)?; + } + if !model.reference_points().is_empty() { + self.comment("Reference points"); + let mut points: Vec<_> = model.reference_points().iter().collect(); + points.sort_by_key(|(_, p)| p.label); + self.line("*NODE"); + for (name, p) in points { + self.comment(&format!("reference point {name}")); + self.numbers_line(&p.label.to_string(), &p.xyz)?; + } + } + Ok(()) + } + + fn elements(&mut self, by_type: &BTreeMap>) { + self.comment("Elements"); + for (t, elements) in by_type { + self.line(&format!("*ELEMENT, TYPE={t}")); + for (e, element) in elements { + let _ = write!(self.out, "{e}"); + for (i, n) in element.nodes.iter().enumerate() { + // A trailing comma continues the data line. + let continues = i >= FIRST_LINE_NODES && (i - FIRST_LINE_NODES).is_multiple_of(CONTINUATION_NODES); + let _ = write!(self.out, "{}{n}", if continues { ",\n" } else { ", " }); + } + self.out.push('\n'); + } + } + } + + fn sets(&mut self, model: &Model, r: &Resolved) -> Result<()> { + if r.sets.is_empty() { + return Ok(()); + } + self.comment("Sets"); + for (name, s) in &r.sets { + if s.is_empty() { + self.comment(&format!("set {name} resolves to nothing on this mesh")); + continue; + } + if !s.nodes.is_empty() { + self.line(&format!("*NSET, NSET={name}")); + self.labels(&s.nodes); + } + if !s.elements.is_empty() { + self.line(&format!("*ELSET, ELSET={name}")); + self.labels(&s.elements); + } + if !s.faces.is_empty() { + self.line(&format!("*SURFACE, TYPE=ELEMENT, NAME={name}")); + for f in &s.faces { + let id = model.element(f.element).and_then(|e| face_id(e.kind, f.face)).ok_or_else(|| Error::Unsupported { + format: FORMAT, + what: format!("set `{name}`: element {} has no face {}", f.element, f.face), + })?; + let _ = writeln!(self.out, "{}, {id}", f.element); + } + } + } + Ok(()) + } + + fn materials(&mut self, model: &Model) -> Result<()> { + if model.materials().is_empty() { + return Ok(()); + } + self.comment("Materials"); + for (name, m) in model.materials() { + self.line(&format!("*MATERIAL, NAME={name}")); + match m.elastic { + Some(Elastic::Isotropic { young, poisson }) => { + self.line("*ELASTIC, TYPE=ISOTROPIC"); + self.numbers_line("", &[young, poisson])?; + } + Some(Elastic::EngineeringConstants { young, poisson, shear }) => { + // E1, E2, E3, ν12, ν13, ν23, G12, G13 on the first line, G23 on the second. + self.line("*ELASTIC, TYPE=ENGINEERING CONSTANTS"); + let [g12, g13, g23] = shear; + self.numbers_line("", &[young[0], young[1], young[2], poisson[0], poisson[1], poisson[2], g12, g13])?; + self.numbers_line("", &[g23])?; + } + None => {} + } + if let Some(rho) = m.density { + self.line("*DENSITY"); + self.numbers_line("", &[rho])?; + } + } + Ok(()) + } + + fn sections(&mut self, model: &Model) -> Result<()> { + if model.sections().is_empty() { + return Ok(()); + } + self.comment("Sections"); + for (name, s) in model.sections() { + self.comment(&format!("section {name}")); + let (set, material) = (&s.set, &s.material); + match s.kind { + SectionKind::Solid {} => { + self.line(&format!("*SOLID SECTION, ELSET={set}, MATERIAL={material}")); + self.line(","); + } + SectionKind::Rod { area } => { + self.line(&format!("*SOLID SECTION, ELSET={set}, MATERIAL={material}")); + self.numbers_line("", &[area])?; + } + SectionKind::Shell { thickness, offset, points } => { + let offset = if offset == 0.0 { String::new() } else { format!(", OFFSET={}", number(offset)?) }; + self.line(&format!("*SHELL SECTION, ELSET={set}, MATERIAL={material}{offset}")); + push_number(&mut self.out, thickness)?; + let _ = writeln!(self.out, ", {points}"); + } + SectionKind::Beam { profile, orientation } => { + // Abaqus profiles: RECT a, b; CIRC r; PIPE r, t; I l, h, b1, b2, t1, t2, t3 + // (l: from the section origin to the bottom; centred here); BOX a, b, t1–t4. + let (shape, dims) = match profile { + BeamProfile::Rectangle { width, height } => ("RECT", vec![width, height]), + BeamProfile::Circle { radius } => ("CIRC", vec![radius]), + BeamProfile::Pipe { radius, thickness } => ("PIPE", vec![radius, thickness]), + BeamProfile::I { height, width, web, flange } => ("I", vec![height / 2.0, height, width, width, flange, flange, web]), + BeamProfile::Box { width, height, thickness } => ("BOX", vec![width, height, thickness, thickness, thickness, thickness]), + }; + self.line(&format!("*BEAM SECTION, SECTION={shape}, ELSET={set}, MATERIAL={material}")); + self.numbers_line("", &dims)?; + self.numbers_line("", &orientation)?; + } + } + } + Ok(()) + } + + fn couplings(&mut self, model: &Model, r: &Resolved) { + if model.couplings().is_empty() { + return; + } + let surfaces = coupling_surfaces(model, r); + self.comment("Couplings: every node of the set weighted equally"); + for (name, c) in model.couplings() { + let (Some(surface), Some(node)) = (surfaces.get(name), r.reference_point_nodes.get(&c.reference_point)) else { continue }; + self.line(&format!("*SURFACE, TYPE=NODE, NAME={surface}")); + self.line(&format!("{}, 1.0", c.set)); + self.line(&format!("*COUPLING, CONSTRAINT NAME={name}, REF NODE={node}, SURFACE={surface}")); + self.line(match c.kind { + CouplingKind::Kinematic => "*KINEMATIC", + CouplingKind::Distributing => "*DISTRIBUTING, WEIGHTING METHOD=UNIFORM", + }); + for (a, b) in dof_runs(&c.dofs.list()) { + self.line(&format!("{a}, {b}")); + } + } + } + + fn boundary_lines<'n>(&mut self, model: &Model, r: &Resolved, names: impl IntoIterator) -> Result<()> { + for name in names { + let Some(bc) = model.boundary_condition(name) else { continue }; + let target = match bc.target() { + Target::Set(s) => s.clone(), + Target::ReferencePoint(p) => match r.reference_point_nodes.get(p) { + Some(n) => n.to_string(), + None => continue, + }, + }; + self.comment(&format!("boundary condition {name}")); + match bc { + BoundaryCondition::Fixed { .. } => self.line(&format!("{target}, ENCASTRE")), + BoundaryCondition::Symmetry { plane, .. } => { + let s = match plane { + Symmetry::X => "XSYMM", + Symmetry::Y => "YSYMM", + Symmetry::Z => "ZSYMM", + }; + self.line(&format!("{target}, {s}")); + } + BoundaryCondition::Displacement { values, .. } => { + for (dof, v) in (1u8..).zip(values) { + if let Some(v) = v { + self.numbers_line(&format!("{target}, {dof}, {dof}"), &[*v])?; + } + } + } + } + } + Ok(()) + } + + /// Boundary conditions of a step against the state it starts from. + fn step_boundaries(&mut self, model: &Model, r: &Resolved, base: &BTreeSet, active: &BTreeSet) -> Result<()> { + if active == base { + self.comment("boundary conditions: as in the base state"); + Ok(()) + } else if active.is_superset(base) { + self.line("*BOUNDARY"); + self.boundary_lines(model, r, active.difference(base)) + } else { + self.line("*BOUNDARY, OP=NEW"); + self.boundary_lines(model, r, active) + } + } + + /// The keyword of one load family with `OP=NEW`, when the step has entries for it or the + /// previous one had some to remove. + fn op_new(&mut self, keyword: &str, now: bool, before: bool) -> bool { + if now || before { + self.line(&format!("{keyword}, OP=NEW")); + } + now + } + + fn step_loads(&mut self, loads: &StepLoads, before: Written) -> Result<()> { + if self.op_new("*CLOAD", !loads.concentrated.is_empty(), before.concentrated) { + for ((node, dof), v) in &loads.concentrated { + self.numbers_line(&format!("{node}, {dof}"), &[*v])?; + } + } + if self.op_new("*DSLOAD", !loads.pressure.is_empty(), before.pressure) { + for (surface, p) in &loads.pressure { + self.numbers_line(&format!("{surface}, P"), &[*p])?; + } + } + if self.op_new("*DLOAD", !loads.element.is_empty(), before.element) { + for (target, kind, v) in &loads.element { + let target = match target { + ElementTarget::Set(s) => s.to_string(), + ElementTarget::Element(e) => e.to_string(), + }; + match kind { + BodyKind::Gravity => { + let g = (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt(); + self.numbers_line(&format!("{target}, GRAV"), &[g, v[0] / g, v[1] / g, v[2] / g])?; + } + BodyKind::Force => { + for (label, c) in ["BX", "BY", "BZ"].into_iter().zip(v) { + if *c != 0.0 { + self.numbers_line(&format!("{target}, {label}"), &[*c])?; + } + } + } + } + } + } + Ok(()) + } + + fn output(&mut self, step: &Step) { + let pick = |names: [(OutputVariable, &'static str); 2]| -> Vec<&'static str> { + names.into_iter().filter(|(v, _)| step.output.contains(v)).map(|(_, n)| n).collect() + }; + let node = pick([(OutputVariable::Displacement, "U"), (OutputVariable::ReactionForce, "RF")]); + let element = pick([(OutputVariable::Stress, "S"), (OutputVariable::Strain, "E")]); + if node.is_empty() && element.is_empty() { + self.line("*OUTPUT, FIELD, OP=NEW, FREQUENCY=0"); + return; + } + self.line("*OUTPUT, FIELD, OP=NEW"); + if !node.is_empty() { + self.line("*NODE OUTPUT"); + self.line(&node.join(", ")); + } + if !element.is_empty() { + self.line("*ELEMENT OUTPUT"); + self.line(&element.join(", ")); + } + } + + fn steps(&mut self, model: &Model, r: &Resolved) -> Result<()> { + let steps = model.steps(); + // The initial state: zero-valued boundary conditions active in every step. + let initial: BTreeSet = if steps.is_empty() { + BTreeSet::new() + } else { + let in_every_step = |n: &String| steps.iter().all(|s| s.boundary_conditions.contains(n)); + model.boundary_conditions().iter().filter(|(n, bc)| holds_at_zero(bc) && in_every_step(n)).map(|(n, _)| n.clone()).collect() + }; + if !initial.is_empty() { + self.comment("Initial state: boundary conditions active in every step"); + self.line("*BOUNDARY"); + self.boundary_lines(model, r, &initial)?; + } + // Gravity on the whole model acts on every element, through the section sets. + let whole_model_gravity = model.loads().values().any(|l| matches!(l, Load::Gravity { set: None, .. })); + let everything: (Vec<&str>, Vec) = if whole_model_gravity { + (model.sections().values().map(|s| s.set.as_str()).collect(), model.elements().keys().copied().collect()) + } else { + Default::default() + }; + let mut base_bcs = initial; + let mut base_loads: BTreeSet = BTreeSet::new(); + let mut written = Written::default(); + for step in steps { + self.comment(&format!("Step {}", step.name)); + match step.procedure { + Procedure::Static {} => { + self.line(&format!("*STEP, NAME={}, NLGEOM=NO", step.name)); + self.line("*STATIC"); + self.line("1.0, 1.0"); + self.step_boundaries(model, r, &base_bcs, &step.boundary_conditions)?; + if step.loads == base_loads { + self.comment("loads: as in the previous general step"); + } else { + let loads = StepLoads::new(model, r, step, &everything); + self.step_loads(&loads, written)?; + written = loads.written(); + } + base_bcs.clone_from(&step.boundary_conditions); + base_loads.clone_from(&step.loads); + } + Procedure::Frequency { modes, min_hz, max_hz } => { + self.line(&format!("*STEP, NAME={}", step.name)); + self.line("*FREQUENCY, EIGENSOLVER=LANCZOS"); + // Modes, then the frequency range; empty trailing fields are left off. + let _ = write!(self.out, "{modes}"); + let range: Vec> = [min_hz, max_hz].into_iter().rev().skip_while(Option::is_none).collect(); + for f in range.into_iter().rev() { + self.out.push_str(", "); + if let Some(f) = f { + push_number(&mut self.out, f)?; + } + } + self.out.push('\n'); + self.step_boundaries(model, r, &base_bcs, &step.boundary_conditions)?; + if !step.loads.is_empty() { + self.comment("loads have no effect on natural frequencies and are not written"); + } + } + } + self.output(step); + self.line("*END STEP"); + } + Ok(()) + } +} diff --git a/crates/decks/src/lib.rs b/crates/decks/src/lib.rs index 3cbeef9..aecb64e 100644 --- a/crates/decks/src/lib.rs +++ b/crates/decks/src/lib.rs @@ -2,9 +2,19 @@ //! keyword files. Each format is its own module (`abaqus`, `nastran`, `lsdyna`), written from the //! vendor's public reference manual in our own words; examples and expected files are our own. //! Writers read [`modelcraft_model::Model`] and its `resolve()` result, never another format. +//! +//! Writers produce text; saving it is the caller's business (the engine's `deck.write`), so the +//! crate stays free of the file system. [`FORMATS`] lists every writer by the name `deck.write` +//! takes. #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] #![forbid(unsafe_code)] +pub mod abaqus; + +use std::collections::{BTreeMap, BTreeSet}; + +use modelcraft_model::{Diagnostic, Load, Model, Resolved}; + #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum Error { /// The model has something this format cannot express (named, with the format). @@ -18,3 +28,63 @@ pub enum Error { } pub type Result = std::result::Result; + +/// A deck format. +pub struct Format { + /// The name commands use (`abaqus`, …). + pub name: &'static str, + /// The usual file extension, without the dot. + pub extension: &'static str, + writer: fn(&Model, &Resolved) -> Result, +} + +/// A written deck: its text and the model check's warnings. +#[derive(Debug, Clone, PartialEq)] +pub struct Deck { + pub text: String, + pub warnings: Vec, +} + +impl Format { + /// Resolve the model once and write it in this format. + pub fn write(&self, model: &Model) -> Result { + let resolved = model.resolve()?; + let text = (self.writer)(model, &resolved)?; + Ok(Deck { text, warnings: resolved.warnings }) + } +} + +/// Every deck format, one row per writer. +pub const FORMATS: &[Format] = &[Format { name: "abaqus", extension: "inp", writer: abaqus::write_resolved }]; + +/// The format called `name`. +pub fn format(name: &str) -> Option<&'static Format> { + FORMATS.iter().find(|f| f.name == name) +} + +/// Fail with [`Error::Unsupported`] naming every item, or succeed when there are none. Writers +/// collect everything their format cannot express and report it at once. +pub(crate) fn unsupported(format: &'static str, items: BTreeSet) -> Result<()> { + if items.is_empty() { Ok(()) } else { Err(Error::Unsupported { format, what: items.into_iter().collect::>().join("; ") }) } +} + +/// The concentrated loads of a set of active loads: the resolved nodal values of every force, +/// traction (degrees of freedom 1–3) and moment (4–6), summed per node and degree of freedom so +/// loads sharing a node add up; zero sums are left out. +pub fn concentrated_loads<'a>(model: &Model, r: &Resolved, active: impl IntoIterator) -> BTreeMap<(u32, u8), f64> { + let mut out: BTreeMap<(u32, u8), f64> = BTreeMap::new(); + for name in active { + let first_dof = match model.load(name) { + Some(Load::Moment { .. }) => 4, + Some(Load::Force { .. } | Load::Traction { .. }) => 1, + _ => continue, + }; + for (node, v) in r.nodal_loads.get(name).map(Vec::as_slice).unwrap_or_default() { + for (dof, c) in (first_dof..).zip(v) { + *out.entry((*node, dof)).or_default() += c; + } + } + } + out.retain(|_, v| *v != 0.0); + out +} diff --git a/crates/decks/tests/abaqus.rs b/crates/decks/tests/abaqus.rs new file mode 100644 index 0000000..1057c2e --- /dev/null +++ b/crates/decks/tests/abaqus.rs @@ -0,0 +1,700 @@ +//! Abaqus writer tests: golden decks of our own small models, a re-read of the written mesh and +//! sets by a minimal test-only reader, element orientation (Jacobian sign) from the written +//! connectivity, number round trips, line limits and everything reported as unsupported. + +use std::collections::{BTreeMap, BTreeSet}; +use std::path::PathBuf; + +use modelcraft_decks::{Error, abaqus}; +use modelcraft_model::{ + BeamProfile, BoundaryCondition, Coupling, CouplingKind, Distribution, Dofs, Elastic, ElementFace, ElementKind, Formulation, Integration, Load, + Material, Model, OutputVariable, Procedure, ReferencePoint, Section, SectionKind, Set, Step, Symmetry, Target, UnitSystem, +}; + +// --------------------------------------------------------------------------------------------- +// Models + +fn steel() -> Material { + Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9) } +} + +fn faces(list: &[(u32, u8)]) -> Set { + Set { faces: list.iter().map(|&(element, face)| ElementFace { element, face }).collect(), ..Set::default() } +} + +fn elements(list: impl IntoIterator) -> Set { + Set { elements: list.into_iter().collect(), ..Set::default() } +} + +fn step(name: &str, procedure: Procedure, bcs: &[&str], loads: &[&str]) -> Step { + let mut s = Step::new(name, procedure); + s.boundary_conditions = bcs.iter().map(|s| s.to_string()).collect(); + s.loads = loads.iter().map(|s| s.to_string()).collect(); + s +} + +/// A 100 × 10 × 10 cantilever of four hex8 elements along X: node (i, j, k) at +/// (25 i, 10 j, 10 k) is labelled 1 + 4 i + 2 j + k. +fn cantilever_mesh() -> Model { + let mut m = Model::default(); + m.set_units(UnitSystem::MmTonneS); + let n = |i: u32, j: u32, k: u32| 1 + 4 * i + 2 * j + k; + for i in 0..=4 { + for j in 0..=1 { + for k in 0..=1 { + m.set_node(n(i, j, k), [25.0 * f64::from(i), 10.0 * f64::from(j), 10.0 * f64::from(k)]).unwrap(); + } + } + } + for i in 0..4 { + let nodes = vec![n(i, 0, 0), n(i + 1, 0, 0), n(i + 1, 1, 0), n(i, 1, 0), n(i, 0, 1), n(i + 1, 0, 1), n(i + 1, 1, 1), n(i, 1, 1)]; + m.set_element(i + 1, ElementKind::Hex8, nodes).unwrap(); + } + m.add_set("all", elements(1..=4)).unwrap(); + // Face 6 of element 1 lies on x = 0, face 4 of element 4 on x = 100. + m.add_set("root", faces(&[(1, 6)])).unwrap(); + m.add_set("tip", faces(&[(4, 4)])).unwrap(); + m.add_material("Steel", steel()).unwrap(); + let formulation = Formulation { incompatible_modes: true, ..Formulation::default() }; + m.add_section("beam", Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation }).unwrap(); + m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("root".into()) }).unwrap(); + m +} + +/// The cantilever with a total end force of 100 down. +fn cantilever() -> Model { + let mut m = cantilever_mesh(); + m.add_load("end", Load::Force { target: Target::Set("tip".into()), vector: [0.0, 0.0, -100.0], distribution: Distribution::Total }).unwrap(); + m.add_step(step("Bend", Procedure::Static {}, &["clamp"], &["end"]), None).unwrap(); + m +} + +/// A 200 × 200 plate of 2 × 2 quad4 shells in z = 0, clamped along x = 0, under pressure on its +/// top face and its own weight, with a symmetry plane along x = 100. +fn plate() -> Model { + let mut m = Model::default(); + m.set_units(UnitSystem::MmTonneS); + let n = |i: u32, j: u32| 1 + i + 3 * j; + for j in 0..=2 { + for i in 0..=2 { + m.set_node(n(i, j), [100.0 * f64::from(i), 100.0 * f64::from(j), 0.0]).unwrap(); + } + } + let mut e = 1; + for j in 0..2 { + for i in 0..2 { + m.set_element(e, ElementKind::Quad4, vec![n(i, j), n(i + 1, j), n(i + 1, j + 1), n(i, j + 1)]).unwrap(); + e += 1; + } + } + m.add_set("skin", elements(1..=4)).unwrap(); + m.add_set("top", faces(&[(1, 1), (2, 1), (3, 1), (4, 1)])).unwrap(); + m.add_set("edge", Set { nodes: vec![1, 4, 7], ..Set::default() }).unwrap(); + m.add_set("mid", Set { nodes: vec![2, 5, 8], ..Set::default() }).unwrap(); + let aluminium = Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9) }; + m.add_material("Aluminium", aluminium).unwrap(); + let kind = SectionKind::Shell { thickness: 2.0, offset: 0.5, points: 5 }; + let formulation = Formulation { integration: Integration::Reduced, ..Formulation::default() }; + m.add_section("plate", Section { set: "skin".into(), material: "Aluminium".into(), kind, formulation }).unwrap(); + m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("edge".into()) }).unwrap(); + m.add_boundary_condition("sym", BoundaryCondition::Symmetry { target: Target::Set("mid".into()), plane: Symmetry::X }).unwrap(); + m.add_load("wind", Load::Pressure { set: "top".into(), magnitude: 0.01 }).unwrap(); + m.add_load("weight", Load::Gravity { acceleration: [0.0, 0.0, -9810.0], set: None }).unwrap(); + m.add_step(step("Load", Procedure::Static {}, &["clamp", "sym"], &["wind", "weight"]), None).unwrap(); + m +} + +/// The cantilever with a reference point coupled to its tip, three static steps that change the +/// active boundary conditions and loads, and a frequency step in between. +fn frequency() -> Model { + let mut m = cantilever_mesh(); + m.add_reference_point("Tip", ReferencePoint { xyz: [110.0, 5.0, 5.0], label: 0 }).unwrap(); + m.add_coupling("tie", Coupling { reference_point: "Tip".into(), set: "tip".into(), kind: CouplingKind::Kinematic, dofs: Dofs::all() }).unwrap(); + let rp = || Target::ReferencePoint("Tip".into()); + m.add_load("twist", Load::Moment { target: rp(), vector: [5000.0, 0.0, 0.0], distribution: Distribution::Total }).unwrap(); + m.add_load("swell", Load::BodyForce { set: "all".into(), vector: [0.0, 1e-3, 0.0] }).unwrap(); + m.add_boundary_condition("lift", BoundaryCondition::Displacement { target: rp(), values: [None, None, Some(0.5), None, None, None] }).unwrap(); + m.add_step(step("Twist", Procedure::Static {}, &["clamp"], &["twist"]), None).unwrap(); + m.add_step(step("Modes", Procedure::Frequency { modes: 6, min_hz: Some(1.0), max_hz: None }, &["clamp", "lift"], &[]), None).unwrap(); + m.add_step(step("Lift", Procedure::Static {}, &["clamp", "lift"], &["twist", "swell"]), None).unwrap(); + let mut release = step("Release", Procedure::Static {}, &["clamp"], &[]); + release.output = [OutputVariable::Displacement].into_iter().collect(); + m.add_step(release, None).unwrap(); + m +} + +// --------------------------------------------------------------------------------------------- +// A minimal reader of what the writer emits (full deck readers are a later milestone) + +#[derive(Default, Debug)] +struct Parsed { + nodes: BTreeMap, + elements: BTreeMap)>, + nsets: BTreeMap>, + elsets: BTreeMap>, + surfaces: BTreeMap>, +} + +fn param(keyword: &str, name: &str) -> Option { + keyword.split(',').skip(1).find_map(|p| { + let (k, v) = p.split_once('=')?; + (k.trim().eq_ignore_ascii_case(name)).then(|| v.trim().to_string()) + }) +} + +fn read(deck: &str) -> Parsed { + let mut p = Parsed::default(); + let mut keyword = String::new(); + let mut pending: Vec = Vec::new(); + for line in deck.lines() { + if line.starts_with("**") { + continue; + } + if line.starts_with('*') { + keyword = line.to_string(); + continue; + } + let head = keyword.split(',').next().unwrap_or_default().trim().to_ascii_uppercase(); + let fields: Vec = line.split(',').map(|f| f.trim().to_string()).collect(); + let label = |s: &str| s.parse::().unwrap(); + match head.as_str() { + "*NODE" => { + let xyz = [1, 2, 3].map(|i| fields[i].parse::().unwrap()); + p.nodes.insert(label(&fields[0]), xyz); + } + "*ELEMENT" => { + pending.extend(fields.iter().filter(|f| !f.is_empty()).cloned()); + if !line.trim_end().ends_with(',') { + let ids: Vec = pending.drain(..).map(|f| label(&f)).collect(); + p.elements.insert(ids[0], (param(&keyword, "TYPE").unwrap(), ids[1..].to_vec())); + } + } + "*NSET" => p.nsets.entry(param(&keyword, "NSET").unwrap()).or_default().extend(fields.iter().map(|f| label(f))), + "*ELSET" => p.elsets.entry(param(&keyword, "ELSET").unwrap()).or_default().extend(fields.iter().map(|f| label(f))), + "*SURFACE" if param(&keyword, "TYPE").as_deref() == Some("ELEMENT") => { + p.surfaces.entry(param(&keyword, "NAME").unwrap()).or_default().push((label(&fields[0]), fields[1].clone())); + } + _ => {} + } + } + p +} + +// --------------------------------------------------------------------------------------------- +// Golden decks + +fn golden(name: &str, deck: &str) { + let path = PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("tests/abaqus").join(format!("{name}.inp")); + if std::env::var_os("MODELCRAFT_BLESS").is_some() { + std::fs::create_dir_all(path.parent().unwrap()).unwrap(); + std::fs::write(&path, deck).unwrap(); + } + let expected = std::fs::read_to_string(&path).unwrap_or_else(|e| panic!("{}: {e} (run with MODELCRAFT_BLESS=1 to write it)", path.display())); + assert!(expected == deck, "{name}.inp differs from the written deck (MODELCRAFT_BLESS=1 rewrites it):\n{deck}"); +} + +fn assert_lines_fit(deck: &str) { + for (i, line) in deck.lines().enumerate() { + assert!(line.len() <= abaqus::MAX_LINE, "line {} has {} characters", i + 1, line.len()); + assert!(line.is_ascii(), "line {} is not ASCII", i + 1); + } +} + +#[test] +fn golden_decks() { + for (name, model) in [("cantilever", cantilever()), ("plate", plate()), ("frequency", frequency())] { + let deck = abaqus::write(&model).unwrap(); + assert_lines_fit(&deck); + golden(name, &deck); + // Deterministic: the same model, and the model after a file round trip, give the same bytes. + assert_eq!(abaqus::write(&model).unwrap(), deck); + let reread: Model = serde_json::from_str(&serde_json::to_string(&model).unwrap()).unwrap(); + assert_eq!(abaqus::write(&reread).unwrap(), deck, "{name}"); + } +} + +#[test] +fn the_end_force_is_shared_by_area_and_keeps_its_total() { + let deck = abaqus::write(&cantilever()).unwrap(); + let cload: Vec<&str> = deck.lines().skip_while(|l| !l.starts_with("*CLOAD")).skip(1).take_while(|l| !l.starts_with('*')).collect(); + // Four equal corners of a flat square face. + assert_eq!(cload, ["17, 3, -25.0", "18, 3, -25.0", "19, 3, -25.0", "20, 3, -25.0"]); + let total: f64 = cload.iter().map(|l| l.rsplit(", ").next().unwrap().parse::().unwrap()).sum(); + assert_eq!(total, -100.0); +} + +#[test] +fn the_written_mesh_and_sets_read_back() { + for model in [cantilever(), plate(), frequency()] { + let deck = abaqus::write(&model).unwrap(); + let p = read(&deck); + let rps: BTreeMap = model.reference_points().values().map(|r| (r.label, r.xyz)).collect(); + let mut nodes: BTreeMap = model.nodes().iter().map(|(l, n)| (*l, n.xyz)).collect(); + nodes.extend(rps); + assert_eq!(p.nodes, nodes); + let connectivity: BTreeMap> = p.elements.iter().map(|(l, (_, n))| (*l, n.clone())).collect(); + assert_eq!(connectivity, model.elements().iter().map(|(l, e)| (*l, e.nodes.clone())).collect()); + let r = model.resolve().unwrap(); + for (name, s) in &r.sets { + assert_eq!(p.nsets.get(name).cloned().unwrap_or_default(), s.nodes, "nset {name}"); + assert_eq!(p.elsets.get(name).cloned().unwrap_or_default(), s.elements, "elset {name}"); + let faces: Vec<(u32, String)> = + s.faces.iter().map(|f| (f.element, abaqus::face_id(model.element(f.element).unwrap().kind, f.face).unwrap().to_string())).collect(); + assert_eq!(p.surfaces.get(name).cloned().unwrap_or_default(), faces, "surface {name}"); + } + } +} + +// --------------------------------------------------------------------------------------------- +// Node order: the element Jacobian of the written connectivity, in Abaqus's own conventions + +/// Natural coordinates of each node, by Abaqus element family. +fn natural(kind: ElementKind) -> Vec<[f64; 3]> { + let hex = [[-1., -1., -1.], [1., -1., -1.], [1., 1., -1.], [-1., 1., -1.], [-1., -1., 1.], [1., -1., 1.], [1., 1., 1.], [-1., 1., 1.]]; + match kind { + ElementKind::Hex8 => hex.to_vec(), + ElementKind::Hex20 => { + let mut v = hex.to_vec(); + // Mid-sides: the bottom face's edges, the top face's, then the vertical ones. + v.extend([[0., -1., -1.], [1., 0., -1.], [0., 1., -1.], [-1., 0., -1.]]); + v.extend([[0., -1., 1.], [1., 0., 1.], [0., 1., 1.], [-1., 0., 1.]]); + v.extend([[-1., -1., 0.], [1., -1., 0.], [1., 1., 0.], [-1., 1., 0.]]); + v + } + ElementKind::Tet4 | ElementKind::Tet10 => { + let mut v = vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]; + if kind == ElementKind::Tet10 { + v.extend([[0.5, 0., 0.], [0.5, 0.5, 0.], [0., 0.5, 0.], [0., 0., 0.5], [0.5, 0., 0.5], [0., 0.5, 0.5]]); + } + v + } + ElementKind::Wedge6 => vec![[0., 0., -1.], [1., 0., -1.], [0., 1., -1.], [0., 0., 1.], [1., 0., 1.], [0., 1., 1.]], + ElementKind::Quad4 => vec![[-1., -1., 0.], [1., -1., 0.], [1., 1., 0.], [-1., 1., 0.]], + ElementKind::Quad8 => { + vec![[-1., -1., 0.], [1., -1., 0.], [1., 1., 0.], [-1., 1., 0.], [0., -1., 0.], [1., 0., 0.], [0., 1., 0.], [-1., 0., 0.]] + } + ElementKind::Tri3 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]], + _ => Vec::new(), + } +} + +/// Abaqus's shape functions of an element family at natural point `q`. +fn shape(kind: ElementKind, q: [f64; 3]) -> Vec { + let nodes = natural(kind); + let [x, y, z] = q; + match kind { + ElementKind::Hex8 => nodes.iter().map(|p| (1. + x * p[0]) * (1. + y * p[1]) * (1. + z * p[2]) / 8.).collect(), + ElementKind::Hex20 => nodes + .iter() + .map(|p| { + let (a, b, c) = (1. + x * p[0], 1. + y * p[1], 1. + z * p[2]); + if p[0] == 0. { + (1. - x * x) * b * c / 4. + } else if p[1] == 0. { + a * (1. - y * y) * c / 4. + } else if p[2] == 0. { + a * b * (1. - z * z) / 4. + } else { + a * b * c * (x * p[0] + y * p[1] + z * p[2] - 2.) / 8. + } + }) + .collect(), + ElementKind::Tet4 => vec![1. - x - y - z, x, y, z], + ElementKind::Tet10 => { + let l = [1. - x - y - z, x, y, z]; + let mut n: Vec = l.iter().map(|li| li * (2. * li - 1.)).collect(); + for (a, b) in [(0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)] { + n.push(4. * l[a] * l[b]); + } + n + } + ElementKind::Wedge6 => { + let t = [1. - x - y, x, y]; + t.iter().map(|ti| ti * (1. - z) / 2.).chain(t.iter().map(|ti| ti * (1. + z) / 2.)).collect() + } + ElementKind::Quad4 => nodes.iter().map(|p| (1. + x * p[0]) * (1. + y * p[1]) / 4.).collect(), + ElementKind::Quad8 => nodes + .iter() + .map(|p| { + if p[0] == 0. { + (1. - x * x) * (1. + y * p[1]) / 2. + } else if p[1] == 0. { + (1. + x * p[0]) * (1. - y * y) / 2. + } else { + (1. + x * p[0]) * (1. + y * p[1]) * (x * p[0] + y * p[1] - 1.) / 4. + } + }) + .collect(), + ElementKind::Tri3 => vec![1. - x - y, x, y], + _ => Vec::new(), + } +} + +/// Columns ∂X/∂ξ, ∂X/∂η, ∂X/∂ζ of the mapping at `q` (central differences). +fn jacobian(kind: ElementKind, xyz: &[[f64; 3]], q: [f64; 3]) -> [[f64; 3]; 3] { + let map = |q: [f64; 3]| { + let n = shape(kind, q); + let mut p = [0.0; 3]; + for (ni, x) in n.iter().zip(xyz) { + for k in 0..3 { + p[k] += ni * x[k]; + } + } + p + }; + let h = 1e-6; + [0, 1, 2].map(|d| { + let (mut a, mut b) = (q, q); + a[d] += h; + b[d] -= h; + let (pa, pb) = (map(a), map(b)); + [0, 1, 2].map(|k| (pa[k] - pb[k]) / (2. * h)) + }) +} + +fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] +} + +fn dot(a: [f64; 3], b: [f64; 3]) -> f64 { + a[0] * b[0] + a[1] * b[1] + a[2] * b[2] +} + +/// The smallest Jacobian determinant over the nodes and the centroid (solids), or the smallest +/// normal component along `up` (shells). +fn min_jacobian(kind: ElementKind, xyz: &[[f64; 3]], up: [f64; 3]) -> f64 { + let nodes = natural(kind); + let centroid = [0, 1, 2].map(|k| nodes.iter().map(|p| p[k]).sum::() / nodes.len() as f64); + nodes + .iter() + .chain([¢roid]) + .map(|q| { + let [a, b, c] = jacobian(kind, xyz, *q); + if kind.dimension() == 3 { dot(cross(a, b), c) } else { dot(cross(a, b), up) } + }) + .fold(f64::INFINITY, f64::min) +} + +/// A one-element model of `kind` with a slanted, stretched, positive-volume placement of the +/// model's reference element (mid-side nodes at edge midpoints). +fn one_element(kind: ElementKind) -> Model { + let corners: Vec<[f64; 3]> = match kind { + ElementKind::Tri3 => 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.]], + ElementKind::Wedge6 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [0., 1., 1.]], + _ => vec![[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]], + }; + let mut xyz = corners.clone(); + xyz.resize(kind.node_count(), [0.0; 3]); + for e in kind.edges() { + if let [a, b, m] = **e { + xyz[m] = [0, 1, 2].map(|i| (corners[a][i] + corners[b][i]) / 2.0); + } + } + // A shear and stretch with a positive determinant, then a shift. + let place = |p: [f64; 3]| [3. * p[0] + 0.5 * p[1] + 10., 2. * p[1] + 0.3 * p[2] - 4., 1.5 * p[2] + 0.2 * p[0] + 7.]; + let mut m = Model::default(); + for (i, p) in xyz.iter().enumerate() { + m.set_node(i as u32 + 1, place(*p)).unwrap(); + } + m.set_element(1, kind, (1..=kind.node_count() as u32).collect()).unwrap(); + m.add_set("all", elements([1])).unwrap(); + m.add_set("base", Set { nodes: vec![1], ..Set::default() }).unwrap(); + m.add_material("Steel", steel()).unwrap(); + let (section, formulation) = if kind.dimension() == 3 { + (SectionKind::Solid {}, Formulation::default()) + } else { + (SectionKind::Shell { thickness: 1.0, offset: 0.0, points: 5 }, Formulation { integration: Integration::Reduced, ..Formulation::default() }) + }; + m.add_section("s", Section { set: "all".into(), material: "Steel".into(), kind: section, formulation }).unwrap(); + m.add_boundary_condition("fix", BoundaryCondition::Fixed { target: Target::Set("base".into()) }).unwrap(); + m.add_step(step("Modes", Procedure::Frequency { modes: 1, min_hz: None, max_hz: None }, &["fix"], &[]), None).unwrap(); + m +} + +#[test] +fn written_connectivity_has_a_positive_jacobian_in_abaqus_conventions() { + let kinds = [ElementKind::Tet4, ElementKind::Tet10, ElementKind::Wedge6, ElementKind::Hex8, ElementKind::Hex20]; + let shells = [ElementKind::Tri3, ElementKind::Quad4, ElementKind::Quad8]; + for kind in kinds.into_iter().chain(shells) { + let model = one_element(kind); + let p = read(&abaqus::write(&model).unwrap()); + let (_, nodes) = &p.elements[&1]; + let xyz: Vec<[f64; 3]> = nodes.iter().map(|n| p.nodes[n]).collect(); + // A shell's positive normal is the side its face 1 (SPOS) looks to: the placed +Z here. + let up = cross([3., 0., 0.2], [0.5, 2., 0.]); + let j = min_jacobian(kind, &xyz, up); + assert!(j > 0.0, "{kind:?}: smallest Jacobian {j}"); + + // The test has teeth: a mirrored linear element shows up (and a misplaced quadratic + // mid-side block, below). + if kind.node_count() == natural(kind).len() + && matches!(kind, ElementKind::Tet4 | ElementKind::Wedge6 | ElementKind::Hex8 | ElementKind::Tri3 | ElementKind::Quad4) + { + let mut mirrored = xyz.clone(); + match kind { + ElementKind::Tet4 | ElementKind::Tri3 => mirrored.swap(1, 2), + ElementKind::Wedge6 => { + mirrored.swap(1, 2); + mirrored.swap(4, 5); + } + ElementKind::Quad4 => mirrored.swap(1, 3), + _ => { + mirrored.swap(1, 3); + mirrored.swap(5, 7); + } + } + assert!(min_jacobian(kind, &mirrored, up) < 0.0, "{kind:?}: mirrored"); + } + } + // Hex20 in the order with the vertical edges' mid-sides before the top face's (Nastran's). + let model = one_element(ElementKind::Hex20); + let p = read(&abaqus::write(&model).unwrap()); + let mut xyz: Vec<[f64; 3]> = p.elements[&1].1.iter().map(|n| p.nodes[n]).collect(); + let (top, vertical) = (xyz[12..16].to_vec(), xyz[16..20].to_vec()); + xyz[12..16].copy_from_slice(&vertical); + xyz[16..20].copy_from_slice(&top); + assert!(min_jacobian(ElementKind::Hex20, &xyz, [0.; 3]) < 0.0, "a misplaced mid-side block is caught"); +} + +#[test] +fn model_faces_are_abaqus_faces() { + // Each face's node positions (1-based, Abaqus element library, our transcription). + let table: [(ElementKind, &[&[usize]]); 3] = [ + (ElementKind::Hex8, &[&[1, 2, 3, 4], &[5, 8, 7, 6], &[1, 5, 6, 2], &[2, 6, 7, 3], &[3, 7, 8, 4], &[4, 8, 5, 1]]), + (ElementKind::Tet4, &[&[1, 2, 3], &[1, 4, 2], &[2, 4, 3], &[3, 4, 1]]), + (ElementKind::Wedge6, &[&[1, 2, 3], &[4, 6, 5], &[1, 4, 5, 2], &[2, 5, 6, 3], &[3, 6, 4, 1]]), + ]; + for (kind, faces) in table { + assert_eq!(kind.faces().len(), faces.len()); + for (i, (ours, theirs)) in kind.faces().iter().zip(faces).enumerate() { + let ours: BTreeSet = ours.iter().map(|p| p + 1).collect(); + let theirs: BTreeSet = theirs.iter().copied().collect(); + assert_eq!(ours, theirs, "{kind:?} face {}", i + 1); + assert_eq!(abaqus::face_id(kind, i as u8 + 1), Some(format!("S{}", i + 1).as_str())); + } + assert_eq!(abaqus::face_id(kind, 0), None); + assert_eq!(abaqus::face_id(kind, faces.len() as u8 + 1), None); + } + assert_eq!(abaqus::face_id(ElementKind::Quad4, 1), Some("SPOS")); + assert_eq!(abaqus::face_id(ElementKind::Tri3, 2), Some("SNEG")); + assert_eq!(abaqus::face_id(ElementKind::Beam2, 1), None); +} + +// --------------------------------------------------------------------------------------------- +// Numbers and limits + +#[test] +fn numbers_read_back_exactly() { + for (v, s) in [(0.1, "0.1"), (210_000.0, "210000.0"), (7.85e-9, "7.85E-9"), (1e300, "1.0E300"), (-0.5, "-0.5"), (1e16, "1.0E16")] { + assert_eq!(abaqus::number(v).unwrap(), s); + } + let mut x: u64 = 0x9E37_79B9_7F4A_7C15; + for _ in 0..100_000 { + x ^= x << 13; + x ^= x >> 7; + x ^= x << 17; + let v = f64::from_bits(x); + if !v.is_finite() { + assert!(abaqus::number(v).is_err()); + continue; + } + let s = abaqus::number(v).unwrap(); + assert_eq!(s.parse::().unwrap().to_bits(), v.to_bits(), "{v:e} → {s}"); + assert!(s.len() <= 24, "{s}"); + assert!(!s.contains('e'), "{s}"); + } + for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + assert!(matches!(abaqus::number(bad), Err(Error::Unsupported { .. }))); + } +} + +#[test] +fn largest_labels_longest_numbers_and_names_fit_a_line() { + let mut m = Model::default(); + let base = modelcraft_model::MAX_LABEL - 30; + let ugly = -1.234_567_890_123_456_7e-300; + let mut xyz = one_element(ElementKind::Hex20).nodes().values().map(|n| n.xyz).collect::>(); + xyz[0] = [ugly, ugly, ugly]; + for (i, p) in xyz.iter().enumerate() { + m.set_node(base + i as u32, *p).unwrap(); + } + m.set_element(modelcraft_model::MAX_LABEL, ElementKind::Hex20, (0..20).map(|i| base + i).collect()).unwrap(); + let long = format!("A{}", "x".repeat(modelcraft_model::MAX_NAME_LEN - 1)); + m.add_set(&long, Set { elements: vec![modelcraft_model::MAX_LABEL], nodes: (0..20).map(|i| base + i).collect(), ..Set::default() }).unwrap(); + m.add_material(&long, Material { elastic: Some(Elastic::Isotropic { young: ugly.abs(), poisson: 0.3 }), density: Some(1.0) }).unwrap(); + m.add_section(&long, Section { set: long.clone(), material: long.clone(), kind: SectionKind::Solid {}, formulation: Formulation::default() }) + .unwrap(); + m.add_reference_point(&long, ReferencePoint { xyz: [ugly; 3], label: 0 }).unwrap(); + m.add_coupling(&long, Coupling { reference_point: long.clone(), set: long.clone(), kind: CouplingKind::Kinematic, dofs: Dofs::all() }).unwrap(); + m.add_boundary_condition(&long, BoundaryCondition::Fixed { target: Target::Set(long.clone()) }).unwrap(); + m.add_step(step(&long, Procedure::Frequency { modes: 100_000, min_hz: Some(ugly.abs()), max_hz: Some(1e300) }, &[&long], &[]), None).unwrap(); + let deck = abaqus::write(&m).unwrap(); + assert_lines_fit(&deck); + let p = read(&deck); + assert_eq!(p.elements[&modelcraft_model::MAX_LABEL].1.len(), 20, "a hex20 continues onto a second line"); + assert_eq!(p.nodes[&base], [ugly; 3]); +} + +// --------------------------------------------------------------------------------------------- +// Step activity, load summing + +fn block<'a>(deck: &'a str, step: &str) -> Vec<&'a str> { + let start = format!("*STEP, NAME={step}"); + deck.lines().skip_while(|l| !l.starts_with(&start)).take_while(|l| *l != "*END STEP").filter(|l| !l.starts_with("**")).collect() +} + +#[test] +fn steps_write_the_activity_table() { + let deck = abaqus::write(&frequency()).unwrap(); + // `clamp` holds at zero in every step: the initial state. + let initial: Vec<&str> = deck.lines().skip_while(|l| !l.starts_with("** Initial state")).skip(1).take(3).collect(); + assert_eq!(initial, ["*BOUNDARY", "** boundary condition clamp", "root, ENCASTRE"]); + let twist = block(&deck, "Twist"); + assert!(twist.contains(&"*CLOAD, OP=NEW") && twist.contains(&"99999999, 4, 5000.0"), "{twist:?}"); + assert!(!twist.iter().any(|l| l.starts_with("*BOUNDARY")), "{twist:?}"); + // The frequency step adds `lift` to its base state and writes no loads. + let modes = block(&deck, "Modes"); + assert!(modes.contains(&"*BOUNDARY") && modes.contains(&"99999999, 3, 3, 0.5"), "{modes:?}"); + assert!(!modes.iter().any(|l| l.contains("LOAD")), "{modes:?}"); + assert!(modes.contains(&"6, 1.0"), "{modes:?}"); + // `Lift` adds `lift` again: the frequency step's change did not carry on. + let lift = block(&deck, "Lift"); + assert!(lift.contains(&"*BOUNDARY") && lift.contains(&"99999999, 3, 3, 0.5"), "{lift:?}"); + assert!(lift.contains(&"*DLOAD, OP=NEW") && lift.contains(&"all, BY, 0.001"), "{lift:?}"); + // `Release` drops `lift` and every load. + let release = block(&deck, "Release"); + assert!(release.contains(&"*BOUNDARY, OP=NEW") && release.contains(&"root, ENCASTRE"), "{release:?}"); + for k in ["*CLOAD, OP=NEW", "*DLOAD, OP=NEW"] { + let at = release.iter().position(|l| *l == k).unwrap_or_else(|| panic!("{k} in {release:?}")); + assert!(release.get(at + 1).is_none_or(|l| l.starts_with('*')), "{k} restates nothing: {release:?}"); + } + assert!(!release.contains(&"*DSLOAD, OP=NEW"), "no pressure before, none now"); + assert!(release.contains(&"U") && !release.contains(&"*ELEMENT OUTPUT"), "{release:?}"); +} + +#[test] +fn loads_sharing_nodes_elements_or_surfaces_add_up() { + let mut m = plate(); + m.add_load("gust", Load::Pressure { set: "top".into(), magnitude: 0.02 }).unwrap(); + m.add_set("half", elements([1, 2])).unwrap(); + m.add_load("extra", Load::Gravity { acceleration: [0.0, 0.0, -9810.0], set: Some("half".into()) }).unwrap(); + m.add_load("p1", Load::Force { target: Target::Set("mid".into()), vector: [3.0, 0.0, 0.0], distribution: Distribution::PerNode }).unwrap(); + m.add_load("p2", Load::Force { target: Target::Set("mid".into()), vector: [1.0, 0.0, 0.0], distribution: Distribution::PerNode }).unwrap(); + for l in ["gust", "extra", "p1", "p2"] { + m.set_load_active("Load", l, true).unwrap(); + } + let deck = abaqus::write(&m).unwrap(); + let lines = block(&deck, "Load"); + assert!(lines.contains(&"top, P, 0.03"), "{lines:?}"); + assert!(lines.contains(&"5, 1, 4.0"), "{lines:?}"); + // Two gravity loads on elements 1 and 2: summed per element. + assert!(lines.contains(&"1, GRAV, 19620.0, 0.0, 0.0, -1.0"), "{lines:?}"); + assert!(lines.contains(&"3, GRAV, 9810.0, 0.0, 0.0, -1.0"), "{lines:?}"); + assert!(!lines.iter().any(|l| l.starts_with("skin, GRAV")), "{lines:?}"); +} + +// --------------------------------------------------------------------------------------------- +// What the format cannot express, and models that are not ready + +fn unsupported(m: &Model) -> String { + match abaqus::write(m) { + Err(Error::Unsupported { format: "Abaqus", what }) => what, + other => panic!("expected unsupported, got {other:?}"), + } +} + +fn with_section(kind: ElementKind, section: SectionKind, integration: Integration) -> Model { + let mut m = one_element(kind); + let formulation = Formulation { integration, ..Formulation::default() }; + m.replace_section("s", Section { set: "all".into(), material: "Steel".into(), kind: section, formulation }).unwrap(); + m +} + +#[test] +fn what_abaqus_cannot_express_is_named() { + let shell = SectionKind::Shell { thickness: 1.0, offset: 0.0, points: 5 }; + let mut tri6 = Model::default(); + for (i, p) in [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0.5, 0., 0.], [0.5, 0.5, 0.], [0., 0.5, 0.]].iter().enumerate() { + tri6.set_node(i as u32 + 1, *p).unwrap(); + } + tri6.set_element(1, ElementKind::Tri6, (1..=6).collect()).unwrap(); + tri6.add_set("all", elements([1])).unwrap(); + tri6.add_material("Steel", steel()).unwrap(); + tri6.add_section("s", Section { set: "all".into(), material: "Steel".into(), kind: shell, formulation: Formulation::default() }).unwrap(); + tri6.add_step(step("Modes", Procedure::Frequency { modes: 1, min_hz: None, max_hz: None }, &[], &[]), None).unwrap(); + assert!(unsupported(&tri6).contains("tri6 shells")); + + assert!(unsupported(&with_section(ElementKind::Quad8, shell, Integration::Full)).contains("full integration of quad8")); + assert!(unsupported(&with_section(ElementKind::Tet10, SectionKind::Solid {}, Integration::Reduced)).contains("reduced integration of tet10")); + assert!(unsupported(&with_section(ElementKind::Wedge6, SectionKind::Solid {}, Integration::Reduced)).contains("wedge6")); + // A linear tetrahedron integrates at one point either way. + assert!(abaqus::write(&with_section(ElementKind::Tet4, SectionKind::Solid {}, Integration::Reduced)).unwrap().contains("TYPE=C3D4\n")); + + let mut case = cantilever(); + case.add_set("ROOT", Set { nodes: vec![1], ..Set::default() }).unwrap(); + let what = unsupported(&case); + assert!(what.contains("sets `ROOT` and `root`") || what.contains("sets `root` and `ROOT`"), "{what}"); + + let mut ortho = plate(); + let carbon = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] }; + ortho.replace_material("Aluminium", Material { elastic: Some(carbon), density: Some(1.6e-9) }).unwrap(); + assert!(unsupported(&ortho).contains("orthotropic material `Aluminium` in a shell section")); + + let mut loose = frequency(); + let dofs = Dofs::try_from(vec![1, 2, 4, 5, 6]).unwrap(); + loose.replace_coupling("tie", Coupling { reference_point: "Tip".into(), set: "tip".into(), kind: CouplingKind::Distributing, dofs }).unwrap(); + assert!(unsupported(&loose).contains("coupling `tie`")); + + // Several problems are named together. + let mut both = with_section(ElementKind::Quad8, shell, Integration::Full); + both.add_set("ALL", Set { nodes: vec![1], ..Set::default() }).unwrap(); + assert!(abaqus::write(&both).unwrap_err().to_string().contains("full integration")); +} + +#[test] +fn orthotropic_solids_beams_and_rods_are_written() { + let mut m = cantilever(); + let carbon = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] }; + m.replace_material("Steel", Material { elastic: Some(carbon), density: None }).unwrap(); + let deck = abaqus::write(&m).unwrap(); + assert!(deck.contains("*ELASTIC, TYPE=ENGINEERING CONSTANTS\n140000.0, 10000.0, 10000.0, 0.3, 0.3, 0.4, 5000.0, 5000.0\n3800.0\n"), "{deck}"); + + // A frame of one beam and one rod. + let mut f = Model::default(); + for (i, p) in [[0., 0., 0.], [1000., 0., 0.], [1000., 0., 500.]].iter().enumerate() { + f.set_node(i as u32 + 1, *p).unwrap(); + } + f.set_element(1, ElementKind::Beam2, vec![1, 2]).unwrap(); + f.set_element(2, ElementKind::Rod2, vec![2, 3]).unwrap(); + f.add_set("beam", elements([1])).unwrap(); + f.add_set("rod", elements([2])).unwrap(); + f.add_set("base", Set { nodes: vec![1, 3], ..Set::default() }).unwrap(); + f.add_material("Steel", steel()).unwrap(); + let i = BeamProfile::I { height: 100.0, width: 50.0, web: 5.0, flange: 8.0 }; + let beam = SectionKind::Beam { profile: i, orientation: [0.0, 1.0, 0.0] }; + f.add_section("b", Section { set: "beam".into(), material: "Steel".into(), kind: beam, formulation: Formulation::default() }).unwrap(); + let rod = SectionKind::Rod { area: 50.0 }; + f.add_section("r", Section { set: "rod".into(), material: "Steel".into(), kind: rod, formulation: Formulation::default() }).unwrap(); + f.add_boundary_condition("pin", BoundaryCondition::Fixed { target: Target::Set("base".into()) }).unwrap(); + f.add_step(step("Modes", Procedure::Frequency { modes: 3, min_hz: None, max_hz: Some(500.0) }, &["pin"], &[]), None).unwrap(); + let deck = abaqus::write(&f).unwrap(); + for expected in [ + "*ELEMENT, TYPE=B31\n1, 1, 2\n", + "*ELEMENT, TYPE=T3D2\n2, 2, 3\n", + "*BEAM SECTION, SECTION=I, ELSET=beam, MATERIAL=Steel\n50.0, 100.0, 50.0, 50.0, 8.0, 8.0, 5.0\n0.0, 1.0, 0.0\n", + "*SOLID SECTION, ELSET=rod, MATERIAL=Steel\n50.0\n", + "3, , 500.0\n", + ] { + assert!(deck.contains(expected), "{expected}\n{deck}"); + } +} + +#[test] +fn a_model_that_is_not_ready_is_refused() { + let mut m = cantilever(); + m.remove_section("beam").unwrap(); + assert!(matches!(abaqus::write(&m), Err(Error::Model(_)))); + assert!(matches!(abaqus::write(&Model::default()), Err(Error::Model(_)))); +} diff --git a/crates/decks/tests/abaqus/cantilever.inp b/crates/decks/tests/abaqus/cantilever.inp new file mode 100644 index 0000000..2cff3f8 --- /dev/null +++ b/crates/decks/tests/abaqus/cantilever.inp @@ -0,0 +1,75 @@ +*HEADING +ModelCraft model, units mm-t-s-N-MPa +** Written by ModelCraft. The numbers are in one consistent unit system and are not converted. +** Nodes +*NODE +1, 0.0, 0.0, 0.0 +2, 0.0, 0.0, 10.0 +3, 0.0, 10.0, 0.0 +4, 0.0, 10.0, 10.0 +5, 25.0, 0.0, 0.0 +6, 25.0, 0.0, 10.0 +7, 25.0, 10.0, 0.0 +8, 25.0, 10.0, 10.0 +9, 50.0, 0.0, 0.0 +10, 50.0, 0.0, 10.0 +11, 50.0, 10.0, 0.0 +12, 50.0, 10.0, 10.0 +13, 75.0, 0.0, 0.0 +14, 75.0, 0.0, 10.0 +15, 75.0, 10.0, 0.0 +16, 75.0, 10.0, 10.0 +17, 100.0, 0.0, 0.0 +18, 100.0, 0.0, 10.0 +19, 100.0, 10.0, 0.0 +20, 100.0, 10.0, 10.0 +** Elements +*ELEMENT, TYPE=C3D8I +1, 1, 5, 7, 3, 2, 6, 8, 4 +2, 5, 9, 11, 7, 6, 10, 12, 8 +3, 9, 13, 15, 11, 10, 14, 16, 12 +4, 13, 17, 19, 15, 14, 18, 20, 16 +** Sets +*NSET, NSET=all +1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 +17, 18, 19, 20 +*ELSET, ELSET=all +1, 2, 3, 4 +*NSET, NSET=root +1, 2, 3, 4 +*SURFACE, TYPE=ELEMENT, NAME=root +1, S6 +*NSET, NSET=tip +17, 18, 19, 20 +*SURFACE, TYPE=ELEMENT, NAME=tip +4, S4 +** Materials +*MATERIAL, NAME=Steel +*ELASTIC, TYPE=ISOTROPIC +210000.0, 0.3 +*DENSITY +7.85E-9 +** Sections +** section beam +*SOLID SECTION, ELSET=all, MATERIAL=Steel +, +** Initial state: boundary conditions active in every step +*BOUNDARY +** boundary condition clamp +root, ENCASTRE +** Step Bend +*STEP, NAME=Bend, NLGEOM=NO +*STATIC +1.0, 1.0 +** boundary conditions: as in the base state +*CLOAD, OP=NEW +17, 3, -25.0 +18, 3, -25.0 +19, 3, -25.0 +20, 3, -25.0 +*OUTPUT, FIELD, OP=NEW +*NODE OUTPUT +U, RF +*ELEMENT OUTPUT +S, E +*END STEP diff --git a/crates/decks/tests/abaqus/frequency.inp b/crates/decks/tests/abaqus/frequency.inp new file mode 100644 index 0000000..079dea7 --- /dev/null +++ b/crates/decks/tests/abaqus/frequency.inp @@ -0,0 +1,123 @@ +*HEADING +ModelCraft model, units mm-t-s-N-MPa +** Written by ModelCraft. The numbers are in one consistent unit system and are not converted. +** Nodes +*NODE +1, 0.0, 0.0, 0.0 +2, 0.0, 0.0, 10.0 +3, 0.0, 10.0, 0.0 +4, 0.0, 10.0, 10.0 +5, 25.0, 0.0, 0.0 +6, 25.0, 0.0, 10.0 +7, 25.0, 10.0, 0.0 +8, 25.0, 10.0, 10.0 +9, 50.0, 0.0, 0.0 +10, 50.0, 0.0, 10.0 +11, 50.0, 10.0, 0.0 +12, 50.0, 10.0, 10.0 +13, 75.0, 0.0, 0.0 +14, 75.0, 0.0, 10.0 +15, 75.0, 10.0, 0.0 +16, 75.0, 10.0, 10.0 +17, 100.0, 0.0, 0.0 +18, 100.0, 0.0, 10.0 +19, 100.0, 10.0, 0.0 +20, 100.0, 10.0, 10.0 +** Reference points +*NODE +** reference point Tip +99999999, 110.0, 5.0, 5.0 +** Elements +*ELEMENT, TYPE=C3D8I +1, 1, 5, 7, 3, 2, 6, 8, 4 +2, 5, 9, 11, 7, 6, 10, 12, 8 +3, 9, 13, 15, 11, 10, 14, 16, 12 +4, 13, 17, 19, 15, 14, 18, 20, 16 +** Sets +*NSET, NSET=all +1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 +17, 18, 19, 20 +*ELSET, ELSET=all +1, 2, 3, 4 +*NSET, NSET=root +1, 2, 3, 4 +*SURFACE, TYPE=ELEMENT, NAME=root +1, S6 +*NSET, NSET=tip +17, 18, 19, 20 +*SURFACE, TYPE=ELEMENT, NAME=tip +4, S4 +** Materials +*MATERIAL, NAME=Steel +*ELASTIC, TYPE=ISOTROPIC +210000.0, 0.3 +*DENSITY +7.85E-9 +** Sections +** section beam +*SOLID SECTION, ELSET=all, MATERIAL=Steel +, +** Couplings: every node of the set weighted equally +*SURFACE, TYPE=NODE, NAME=tie-nodes +tip, 1.0 +*COUPLING, CONSTRAINT NAME=tie, REF NODE=99999999, SURFACE=tie-nodes +*KINEMATIC +1, 6 +** Initial state: boundary conditions active in every step +*BOUNDARY +** boundary condition clamp +root, ENCASTRE +** Step Twist +*STEP, NAME=Twist, NLGEOM=NO +*STATIC +1.0, 1.0 +** boundary conditions: as in the base state +*CLOAD, OP=NEW +99999999, 4, 5000.0 +*OUTPUT, FIELD, OP=NEW +*NODE OUTPUT +U, RF +*ELEMENT OUTPUT +S, E +*END STEP +** Step Modes +*STEP, NAME=Modes +*FREQUENCY, EIGENSOLVER=LANCZOS +6, 1.0 +*BOUNDARY +** boundary condition lift +99999999, 3, 3, 0.5 +*OUTPUT, FIELD, OP=NEW +*NODE OUTPUT +U +*END STEP +** Step Lift +*STEP, NAME=Lift, NLGEOM=NO +*STATIC +1.0, 1.0 +*BOUNDARY +** boundary condition lift +99999999, 3, 3, 0.5 +*CLOAD, OP=NEW +99999999, 4, 5000.0 +*DLOAD, OP=NEW +all, BY, 0.001 +*OUTPUT, FIELD, OP=NEW +*NODE OUTPUT +U, RF +*ELEMENT OUTPUT +S, E +*END STEP +** Step Release +*STEP, NAME=Release, NLGEOM=NO +*STATIC +1.0, 1.0 +*BOUNDARY, OP=NEW +** boundary condition clamp +root, ENCASTRE +*CLOAD, OP=NEW +*DLOAD, OP=NEW +*OUTPUT, FIELD, OP=NEW +*NODE OUTPUT +U +*END STEP diff --git a/crates/decks/tests/abaqus/plate.inp b/crates/decks/tests/abaqus/plate.inp new file mode 100644 index 0000000..6c8e201 --- /dev/null +++ b/crates/decks/tests/abaqus/plate.inp @@ -0,0 +1,67 @@ +*HEADING +ModelCraft model, units mm-t-s-N-MPa +** Written by ModelCraft. The numbers are in one consistent unit system and are not converted. +** Nodes +*NODE +1, 0.0, 0.0, 0.0 +2, 100.0, 0.0, 0.0 +3, 200.0, 0.0, 0.0 +4, 0.0, 100.0, 0.0 +5, 100.0, 100.0, 0.0 +6, 200.0, 100.0, 0.0 +7, 0.0, 200.0, 0.0 +8, 100.0, 200.0, 0.0 +9, 200.0, 200.0, 0.0 +** Elements +*ELEMENT, TYPE=S4R +1, 1, 2, 5, 4 +2, 2, 3, 6, 5 +3, 4, 5, 8, 7 +4, 5, 6, 9, 8 +** Sets +*NSET, NSET=edge +1, 4, 7 +*NSET, NSET=mid +2, 5, 8 +*NSET, NSET=skin +1, 2, 3, 4, 5, 6, 7, 8, 9 +*ELSET, ELSET=skin +1, 2, 3, 4 +*NSET, NSET=top +1, 2, 3, 4, 5, 6, 7, 8, 9 +*SURFACE, TYPE=ELEMENT, NAME=top +1, SPOS +2, SPOS +3, SPOS +4, SPOS +** Materials +*MATERIAL, NAME=Aluminium +*ELASTIC, TYPE=ISOTROPIC +70000.0, 0.33 +*DENSITY +2.7E-9 +** Sections +** section plate +*SHELL SECTION, ELSET=skin, MATERIAL=Aluminium, OFFSET=0.5 +2.0, 5 +** Initial state: boundary conditions active in every step +*BOUNDARY +** boundary condition clamp +edge, ENCASTRE +** boundary condition sym +mid, XSYMM +** Step Load +*STEP, NAME=Load, NLGEOM=NO +*STATIC +1.0, 1.0 +** boundary conditions: as in the base state +*DSLOAD, OP=NEW +top, P, 0.01 +*DLOAD, OP=NEW +skin, GRAV, 9810.0, 0.0, 0.0, -1.0 +*OUTPUT, FIELD, OP=NEW +*NODE OUTPUT +U, RF +*ELEMENT OUTPUT +S, E +*END STEP diff --git a/crates/engine/Cargo.toml b/crates/engine/Cargo.toml index a8cb291..ad03d96 100644 --- a/crates/engine/Cargo.toml +++ b/crates/engine/Cargo.toml @@ -13,6 +13,7 @@ modelcraft-model = { workspace = true } modelcraft-geometry = { workspace = true } modelcraft-mesh = { workspace = true } modelcraft-io = { workspace = true } +modelcraft-decks = { workspace = true } serde = { workspace = true } serde_json = { workspace = true } thiserror = { workspace = true } diff --git a/crates/engine/src/cmd/deck.rs b/crates/engine/src/cmd/deck.rs new file mode 100644 index 0000000..53d1932 --- /dev/null +++ b/crates/engine/src/cmd/deck.rs @@ -0,0 +1,65 @@ +//! Solver decks: the model written in a solver's input format. + +use std::path::Path; + +use serde_json::{Value, json}; + +use super::req; +use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; + +pub const WRITE: CommandSpec = CommandSpec { + id: "deck.write", + label: "Write Solver Deck", + params: r#"{"format": "abaqus", "path": "part.inp"} — writes the model as a solver deck (formats: modelcraft_decks::FORMATS); fails naming every error of the model check and everything the format cannot express"#, + kind: CommandKind::Session, + run: write, +}; + +fn write(s: &mut Session, p: &Value) -> Result { + let name: &str = req(p, "format")?; + let format = modelcraft_decks::format(name).ok_or_else(|| { + let shown: String = name.chars().take(32).collect(); + let known: Vec<&str> = modelcraft_decks::FORMATS.iter().map(|f| f.name).collect(); + EngineError::Params(format!("unknown deck format `{shown}`; known: {}", known.join(", "))) + })?; + let path = super::path_param(p, "path")?; + let deck = format.write(&s.model)?; + modelcraft_io::fs::write_atomic(Path::new(path), deck.text.as_bytes(), false)?; + Ok(json!({"format": format.name, "path": path, "bytes": deck.text.len(), "warnings": deck.warnings})) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn writes_the_defined_model_and_refuses_what_it_cannot() { + let dir = std::env::temp_dir().join(format!("modelcraft-deck-write-{}", std::process::id())); + std::fs::create_dir_all(&dir).unwrap(); + let path = dir.join("bar.inp"); + let path_str = path.to_string_lossy().into_owned(); + + let mut s = Session::default(); + crate::cmd::define::tests::defined_bar(&mut s); + let journal = s.journal.len(); + let out = s.execute("deck.write", &json!({"format": "abaqus", "path": path_str})).unwrap(); + let written = std::fs::read_to_string(&path).unwrap(); + assert_eq!(written, modelcraft_decks::abaqus::write(&s.model).unwrap()); + assert_eq!(out["bytes"], written.len()); + assert_eq!(out["warnings"].as_array().map(Vec::len), Some(1), "`press` is in no step: {out}"); + assert!(written.contains("*ELEMENT, TYPE=C3D8I"), "{written}"); + assert_eq!(s.journal.len(), journal, "writing a deck does not change the model"); + assert!(std::fs::read_dir(&dir).unwrap().count() == 1, "no partial file is left behind"); + + // Unknown formats, missing parameters and models that are not ready fail without writing. + let other = dir.join("other.inp").to_string_lossy().into_owned(); + assert!(s.execute("deck.write", &json!({"format": "nope", "path": other})).unwrap_err().to_string().contains("known: abaqus")); + assert!(s.execute("deck.write", &json!({"path": other})).is_err()); + assert!(s.execute("deck.write", &json!({"format": "abaqus"})).is_err()); + let mut empty = Session::default(); + assert!(matches!(empty.execute("deck.write", &json!({"format": "abaqus", "path": other})), Err(EngineError::Deck(_)))); + assert!(!dir.join("other.inp").exists()); + assert!(matches!(s.execute("deck.write", &json!({"format": "abaqus", "path": dir.to_string_lossy()})), Err(EngineError::File(_)))); + let _ = std::fs::remove_dir_all(&dir); + } +} diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index 5e151cb..3ec2311 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -1,5 +1,6 @@ //! The command registry. Each area has its own module; [`ALL`] lists every command. +pub(crate) mod deck; pub(crate) mod define; mod edit; pub(crate) mod file; @@ -66,6 +67,7 @@ pub static ALL: &[CommandSpec] = &[ file::RECOVERY_DISCARD, edit::UNDO, edit::REDO, + deck::WRITE, ]; /// A required label parameter (1..=u32::MAX; the model checks its own range). diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 7e320cc..a937afa 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -40,6 +40,8 @@ pub enum EngineError { Mesh(#[from] modelcraft_mesh::Error), #[error(transparent)] File(#[from] modelcraft_io::Error), + #[error(transparent)] + Deck(#[from] modelcraft_decks::Error), /// A command panicked; the model is unchanged. #[error("internal error in {command}: {message}")] Internal { command: String, message: String }, @@ -229,6 +231,14 @@ mod tests { json!({"path": "/dev/zero"}), json!({"path": "/no/such/design.solidcraft"}), json!({"path": "x".repeat(100_000)}), + json!({"format": "abaqus"}), + json!({"format": "abaqus", "path": ""}), + json!({"format": "abaqus", "path": "/"}), + json!({"format": "abaqus", "path": "/no/such/dir/part.inp"}), + json!({"format": "abaqus", "path": format!("/no/such/{}", "x".repeat(100_000))}), + json!({"format": "abaqus", "path": 7}), + json!({"format": "x".repeat(100_000), "path": "/no/such/part.inp"}), + json!({"format": null, "path": "/no/such/part.inp"}), json!({"path": 7}), json!({"body": "Body1", "tolerance": 0.05}), json!({"body": null, "tolerance": "fine"}), diff --git a/docs/architecture.md b/docs/architecture.md index 6b0deb2..693bf77 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -22,7 +22,7 @@ crates in strictly lower layers. | L1 | `modelcraft-geometry` | The only door to SolidCraft: opens `.solidcraft` designs (rebuilt from their timeline, bodies with persistent face names), B-rep bodies, measures, tessellation into a face-tagged `SurfaceMesh`; face-tagged 3MF files read into the same type without the kernel (docs/formats/tagged-mesh.md). No SolidCraft type in its API. | started | | L1 | `modelcraft-mesh` | Surface preparation (weld, watertightness, feature edges), robust predicates, Delaunay tetrahedralization and refinement, Voronoi-based optimization, hex meshing, quality metrics. | started | | L1 | `modelcraft-solver` | Element formulations, assembly, sparse linear static and modal solution. | planned | -| L1 | `modelcraft-decks` | Abaqus input, Nastran bulk data and LS-DYNA keyword readers and writers, to and from `model`. | planned | +| L1 | `modelcraft-decks` | Abaqus input, Nastran bulk data and LS-DYNA keyword readers and writers, to and from `model`. Text out, no file system; `deck.write` saves it. | started | | L1 | `modelcraft-render` | Display data without a GPU dependency: the visible faces, element edges and nodes of a model, the orbit camera, a CPU rasterizer for headless snapshots; contour maps later. | started | | L2 | `modelcraft-doc` | The document: geometry, mesh, model, results, history. | planned | | L3 | `modelcraft-io` | `.modelcraft` files (docs/formats/modelcraft.md): text model definition, binary mesh, geometry references with content hashes, the journal; atomic saves; autosave storage for crash recovery. Bytes in, bytes out; the file system only in `fs` and `recovery`. Result files later. | started | diff --git a/docs/parity.md b/docs/parity.md index 4a9da46..355ac30 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -56,7 +56,7 @@ products are run, no documentation text is copied). Status: **done**, **started* | 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 | -| Decks | Abaqus input writer | M3 | | planned | +| Decks | Abaqus input writer (`deck.write`, flat deck; reports what Abaqus cannot express) | M3 | | done | | Decks | Nastran bulk data writer (SOL 101, 103) | M3 | | planned | | Decks | LS-DYNA keyword writer (implicit; explicit in M8) | M3 | | planned | | Decks | Readers keeping unknown keywords | M9 | | planned | @@ -76,3 +76,8 @@ products are run, no documentation text is copied). Status: **done**, **started* as degenerate quads. Every mapping is tested by the sign of the element Jacobian after writing. - Each format's field widths, continuation rules, number forms and label limits are respected by the writer and checked by golden-file tests on our own examples. +- Abaqus: the model's node order and face numbers are Abaqus's, so connectivity passes through; + element types come from the element kind, the section and its formulation (table in + `crates/decks/src/abaqus.rs`). The step table is bridged to Abaqus's propagation: zero-valued + boundary conditions active in every step go to the initial state, and steps restate their + boundary conditions and loads with `OP=NEW` when the active set changes.