OneCAD-Simulation/crates/decks/src/abaqus.rs

//! 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<String> {
    write_resolved(model, &model.resolve()?)
}

/// [`write`] for a model already resolved.
pub fn write_resolved(model: &Model, r: &Resolved) -> Result<String> {
    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<String, String> {
    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 (&section.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<BTreeMap<String, Vec<(u32, &'m Element)>>> {
    let mut by_type: BTreeMap<String, Vec<(u32, &Element)>> = BTreeMap::new();
    let mut bad = BTreeSet::new();
    let mut count = 0;
    for (name, section) in model.sections() {
        let mut by_kind: BTreeMap<ElementKind, Vec<(u32, &Element)>> = BTreeMap::new();
        for &e in r.sets.get(&section.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<String, &str> = 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(&section.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<String> {
    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<String, String> {
    let mut taken: BTreeSet<String> = 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<u32> = 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<u32, [f64; 3]> = 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<u32>)) -> 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<String, Vec<(u32, &Element)>>) {
        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<Item = &'n String>) -> 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<String>, active: &BTreeSet<String>) -> 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<String> = 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<u32>) = 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<String> = 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<Option<f64>> = [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(())
    }
}