Nastran bulk data writer (SOL 101/103) and deck.write steps (M3) (#24)

* Nastran bulk data writer: fixed-format entries, exact large-field reals, SOL 101/103 decks

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* Nastran writer tests: golden decks, a re-reader of our output, node order against Nastran's conventions

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* docs: Nastran writer in the roadmap, parity table and deck notes

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* Nastran writer: simplify planning (subcases, one property table), stream cards into the deck, write through the VFS

ElementKind::corner_count/is_quadratic and Dofs::bits/from_bits join the model.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* Nastran joins the shared deck formats; deck.write chooses steps for every format

nastran::write_resolved is the FORMATS row; Format::write_steps and with_steps write a subset of
steps (a Nastran run is static or frequency steps); problems go through the shared unsupported().

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

---------

Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com>
Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

ATTRIBUTION.md

@@ -14,6 +14,7 @@ is committed and listed here with its exports.
 
 | File | Author | Source | Licence |
 |---|---|---|---|
+| crates/decks/tests/golden/nastran/*.bdf | ModelCraft contributors | Written by our Nastran writer from our own test models (crates/decks/tests/nastran.rs) | MIT OR Apache-2.0 |
 
 ## Adapted code
 

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 (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 |
+| M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input (done: `deck.write`), Nastran SOL 101/103 (done: `deck.write`) and LS-DYNA implicit; 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 |

crates/decks/src/lib.rs

@@ -10,6 +10,7 @@
 #![forbid(unsafe_code)]
 
 pub mod abaqus;
+pub mod nastran;
 
 use std::collections::{BTreeMap, BTreeSet};
 
@@ -52,10 +53,40 @@ impl Format {
         let text = (self.writer)(model, &resolved)?;
         Ok(Deck { text, warnings: resolved.warnings })
     }
+
+    /// Write only the named steps (in the model's order), or every step when `steps` is `None`;
+    /// for runs that hold some steps only, such as a Nastran run (static or frequency steps).
+    pub fn write_steps(&self, model: &Model, steps: Option<&[String]>) -> Result<Deck> {
+        match steps {
+            None => self.write(model),
+            Some(names) => self.write(&with_steps(model, names)?),
+        }
+    }
+}
+
+/// The model with only the named steps, in its own order. Definitions the dropped steps used
+/// stay (the check warns about the ones no step uses).
+pub fn with_steps(model: &Model, names: &[String]) -> Result<Model> {
+    if names.is_empty() {
+        return Err(Error::Model(modelcraft_model::Error::Invalid("choose at least one step to write".into())));
+    }
+    if let Some(n) = names.iter().find(|n| model.step(n).is_none()) {
+        let shown: String = n.chars().take(modelcraft_model::MAX_NAME_LEN).collect();
+        return Err(Error::Model(modelcraft_model::Error::NotFound(format!("step `{shown}`"))));
+    }
+    let mut out = model.clone();
+    let dropped: Vec<String> = model.steps().iter().map(|s| s.name.clone()).filter(|s| !names.contains(s)).collect();
+    for name in dropped {
+        out.remove_step(&name)?;
+    }
+    Ok(out)
 }
 
 /// Every deck format, one row per writer.
-pub const FORMATS: &[Format] = &[Format { name: "abaqus", extension: "inp", writer: abaqus::write_resolved }];
+pub const FORMATS: &[Format] = &[
+    Format { name: "abaqus", extension: "inp", writer: abaqus::write_resolved },
+    Format { name: "nastran", extension: "bdf", writer: nastran::write_resolved },
+];
 
 /// The format called `name`.
 pub fn format(name: &str) -> Option<&'static Format> {

crates/decks/src/nastran/card.rs

@@ -0,0 +1,338 @@
+//! Bulk data entries in fixed format, and the numbers in them.
+//!
+//! An entry is a name and a stream of data fields. In small-field format each line holds the name
+//! (or, on a continuation line, blanks) in columns 1–8 and eight data fields of eight columns; in
+//! large-field format the name carries a `*`, the data fields are sixteen columns wide, four per
+//! line, and each continuation line starts with `*`. Continuations are automatic: the next line
+//! continues the entry when its first field is blank (small) or `*` (large), so no continuation
+//! markers are written in field 10. A data stream is the same in both formats: field `k` of an
+//! entry sits at the same place in the stream whichever width carries it.
+//!
+//! Entries that hold only integers and keywords use small fields (labels stop at 99 999 999, eight
+//! digits). Entries that hold a real number use large fields, because eight columns keep at most
+//! seven significant digits. A real is written exactly when its shortest round-trip form fits
+//! sixteen columns (every hand-entered value and coordinates on decimal grids); otherwise it is
+//! rounded to the nearest value that fits, which keeps at least ten significant digits and, for
+//! magnitudes between 1e-5 and 1e9 without a sign, fourteen. Free-field format would not help:
+//! readers round free-field reals to the same widths.
+
+use std::fmt::Write as _;
+
+use crate::{Error, Result};
+
+/// Columns of a small field.
+const SMALL: usize = 8;
+/// Columns of a large field.
+const LARGE: usize = 16;
+
+/// One data field.
+#[derive(Debug, Clone, PartialEq)]
+pub(crate) enum Field {
+    Blank,
+    Int(i64),
+    Real(f64),
+    Text(&'static str),
+}
+
+impl From<u32> for Field {
+    fn from(v: u32) -> Field {
+        Field::Int(v.into())
+    }
+}
+
+impl From<f64> for Field {
+    fn from(v: f64) -> Field {
+        Field::Real(v)
+    }
+}
+
+impl From<&'static str> for Field {
+    fn from(v: &'static str) -> Field {
+        Field::Text(v)
+    }
+}
+
+/// A bulk data entry: a name and its data fields in order (fields 2–9 of the first line, then
+/// fields 2–9 of each continuation).
+#[derive(Debug, Clone, PartialEq)]
+pub(crate) struct Card {
+    name: &'static str,
+    fields: Vec<Field>,
+}
+
+impl Card {
+    pub(crate) fn new(name: &'static str) -> Card {
+        Card { name, fields: Vec::new() }
+    }
+
+    pub(crate) fn int(mut self, v: impl Into<i64>) -> Card {
+        self.fields.push(Field::Int(v.into()));
+        self
+    }
+
+    pub(crate) fn ints<I: Into<i64>>(mut self, v: impl IntoIterator<Item = I>) -> Card {
+        self.fields.extend(v.into_iter().map(|i| Field::Int(i.into())));
+        self
+    }
+
+    pub(crate) fn real(mut self, v: f64) -> Card {
+        self.fields.push(Field::Real(v));
+        self
+    }
+
+    pub(crate) fn reals(mut self, v: impl IntoIterator<Item = f64>) -> Card {
+        self.fields.extend(v.into_iter().map(Field::Real));
+        self
+    }
+
+    pub(crate) fn text(mut self, v: &'static str) -> Card {
+        self.fields.push(Field::Text(v));
+        self
+    }
+
+    pub(crate) fn blank(mut self) -> Card {
+        self.fields.push(Field::Blank);
+        self
+    }
+
+    /// An optional field: blank when absent.
+    pub(crate) fn opt(mut self, v: Option<impl Into<Field>>) -> Card {
+        self.fields.push(v.map_or(Field::Blank, Into::into));
+        self
+    }
+
+    pub(crate) fn name(&self) -> &'static str {
+        self.name
+    }
+
+    /// Blank fields up to the end of the current line, so the next field starts a continuation.
+    pub(crate) fn end_line(mut self) -> Card {
+        while !self.fields.len().is_multiple_of(8) {
+            self.fields.push(Field::Blank);
+        }
+        self
+    }
+
+    /// Append the entry to `out`. Trailing blank fields are dropped, and so are trailing blanks on
+    /// each line.
+    pub(crate) fn write(&self, out: &mut String) -> Result<()> {
+        let large = self.fields.iter().any(|f| matches!(f, Field::Real(_)));
+        let (width, per_line) = if large { (LARGE, 4) } else { (SMALL, 8) };
+        let used = self.fields.iter().rposition(|f| *f != Field::Blank).map_or(0, |i| i + 1);
+        let too_wide = |text: &dyn std::fmt::Display| Error::Limit(format!("{}: `{text}` does not fit a {width}-column field", self.name));
+        let mut lines = self.fields.get(..used).unwrap_or_default().chunks(per_line).peekable();
+        let mut first = true;
+        loop {
+            let start = out.len();
+            if first {
+                out.push_str(self.name);
+            }
+            if large {
+                out.push('*');
+            }
+            while out.len() < start + SMALL {
+                out.push(' ');
+            }
+            for f in lines.next().unwrap_or_default() {
+                match *f {
+                    Field::Blank => out.extend(std::iter::repeat_n(' ', width)),
+                    Field::Int(v) if digits(v) <= width => {
+                        let _ = write!(out, "{v:>width$}");
+                    }
+                    Field::Int(v) => return Err(too_wide(&v)),
+                    Field::Text(t) if t.len() <= width => {
+                        let _ = write!(out, "{t:>width$}");
+                    }
+                    Field::Text(t) => return Err(too_wide(&t)),
+                    Field::Real(v) => {
+                        let (text, _) = real(v, width).ok_or_else(|| Error::Limit(format!("{}: {v} is not a finite number", self.name)))?;
+                        let _ = write!(out, "{text:>width$}");
+                    }
+                }
+            }
+            out.truncate(start + out.get(start..).map_or(0, |l| l.trim_end().len()));
+            out.push('\n');
+            first = false;
+            if lines.peek().is_none() {
+                return Ok(());
+            }
+        }
+    }
+}
+
+/// Characters of an integer written in decimal.
+fn digits(v: i64) -> usize {
+    let sign = usize::from(v < 0);
+    sign + v.unsigned_abs().checked_ilog10().map_or(1, |l| l as usize + 1)
+}
+
+/// A finite real in at most `width` columns, and whether it is exact (reads back as the same
+/// `f64`). `None` for NaN and infinities.
+///
+/// Forms, in order of preference: plain decimal with a leading digit (`0.3`, `210000.`) for
+/// magnitudes from 1e-5 below 1e16; an exponent with `E` (`7.85E-9`); the exponent without `E` (`7.85-9`);
+/// a plain decimal without its leading zero (`.0001`). When none fits the shortest round-trip
+/// digits, the value is rounded to fewer significant digits until one does.
+pub(crate) fn real(v: f64, width: usize) -> Option<(String, bool)> {
+    if !v.is_finite() {
+        return None;
+    }
+    if v == 0.0 {
+        return Some(("0.".to_string(), true));
+    }
+    let (digits, exp) = shortest(v.abs());
+    if let Some(s) = fit(v < 0.0, &digits, exp, width) {
+        return Some((s, true));
+    }
+    for keep in (1..digits.len()).rev() {
+        let rounded: f64 = format!("{:.*e}", keep - 1, v.abs()).parse().ok()?;
+        let (digits, exp) = shortest(rounded);
+        if let Some(s) = fit(v < 0.0, &digits, exp, width) {
+            return Some((s, false));
+        }
+    }
+    None
+}
+
+/// The shortest round-trip significant digits of a positive finite number, and the decimal
+/// exponent of the first digit (`v = d.ddd × 10^exp`).
+fn shortest(v: f64) -> (String, i32) {
+    let s = format!("{v:e}");
+    let (mantissa, exp) = s.split_once('e').unwrap_or((&s, "0"));
+    (mantissa.replace('.', ""), exp.parse().unwrap_or(0))
+}
+
+/// The first form of `±d.ddd × 10^exp` that fits `width` (lengths are worked out before
+/// anything is formatted).
+fn fit(negative: bool, digits: &str, exp: i32, width: usize) -> Option<String> {
+    let sign = if negative { "-" } else { "" };
+    let n = digits.len();
+    let e = exp.unsigned_abs() as usize;
+    let exp_len = e.checked_ilog10().map_or(1, |l| l as usize + 1);
+    // Plain decimal: the integer digits (padded with zeros), or zeros after the point.
+    let plain_len = |leading_zero: bool| sign.len() + if exp >= 0 { n.max(e + 1) + 1 } else { usize::from(leading_zero) + 1 + (e - 1) + n };
+    let plain = |leading_zero: bool| {
+        if exp >= 0 {
+            match digits.split_at_checked(e + 1) {
+                Some((a, b)) => format!("{sign}{a}.{b}"),
+                None => format!("{sign}{digits}{}.", "0".repeat(e + 1 - n)),
+            }
+        } else {
+            let zero = if leading_zero { "0" } else { "" };
+            format!("{sign}{zero}.{}{digits}", "0".repeat(e - 1))
+        }
+    };
+    let scientific = |with_e: bool| {
+        let (first, rest) = digits.split_at(1.min(n));
+        let marker = if with_e { "E" } else { "" };
+        let exp_sign = if exp < 0 { '-' } else { '+' };
+        format!("{sign}{first}.{rest}{marker}{exp_sign}{e}")
+    };
+    let sci_len = sign.len() + n + 1 + 1 + exp_len;
+    if (-5..16).contains(&exp) && plain_len(true) <= width {
+        Some(plain(true))
+    } else if sci_len < width {
+        Some(scientific(true))
+    } else if sci_len <= width {
+        Some(scientific(false))
+    } else if exp < 0 && plain_len(false) <= width {
+        Some(plain(false))
+    } else {
+        None
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// Read a real the way a Nastran reader does: an exponent may come without `E`.
+    fn read(s: &str) -> f64 {
+        let s = s.trim();
+        let split = s.char_indices().skip(1).find(|(i, c)| (*c == '+' || *c == '-') && !s[..*i].ends_with(['E', 'e']));
+        let normal = match split {
+            Some((i, _)) => format!("{}E{}", &s[..i], &s[i..]),
+            None => s.to_string(),
+        };
+        normal.parse().unwrap()
+    }
+
+    #[test]
+    fn reals_take_the_plainest_form_that_fits() {
+        let cases = [
+            (0.0, "0."),
+            (-0.0, "0."),
+            (0.3, "0.3"),
+            (-1.5, "-1.5"),
+            (210000.0, "210000."),
+            (7.85e-9, "7.85E-9"),
+            (1e-5, "0.00001"),
+            (1e20, "1.E+20"),
+            (-2.5e-300, "-2.5E-300"),
+            (123456789.12345, "123456789.12345"),
+            (f64::MAX, "1.7976931348623157E+308"),
+        ];
+        for (v, want) in cases {
+            let (s, exact) = real(v, 32).unwrap();
+            assert_eq!(s, want, "{v}");
+            assert!(exact);
+            assert_eq!(read(&s), v);
+        }
+        assert_eq!(real(f64::NAN, 16), None);
+        assert_eq!(real(f64::INFINITY, 16), None);
+        // Tight widths fall back to the E-less exponent, then to fewer digits.
+        assert_eq!(real(1.25e-12, 7).unwrap(), ("1.25-12".into(), true));
+        assert_eq!(real(1.0 / 3.0, 8).unwrap(), (".3333333".into(), false));
+    }
+
+    #[test]
+    fn sixteen_columns_round_trip_short_values_and_keep_ten_digits_of_any() {
+        let mut seed = 0x9e37_79b9_7f4a_7c15_u64;
+        let mut next = || {
+            seed ^= seed << 13;
+            seed ^= seed >> 7;
+            seed ^= seed << 17;
+            seed
+        };
+        for _ in 0..20_000 {
+            let v = f64::from_bits(next());
+            if !v.is_finite() {
+                continue;
+            }
+            let (s, exact) = real(v, LARGE).unwrap();
+            assert!(s.len() <= LARGE, "{s}");
+            let back = read(&s);
+            if exact {
+                assert_eq!(back.to_bits(), v.to_bits(), "{v} -> {s}");
+            }
+            assert!(((back - v) / v).abs() < 1e-9 || v == 0.0, "{v} -> {s}");
+            // Coordinates and moduli in everyday ranges keep fourteen digits.
+            let w = (next() as f64 / u64::MAX as f64) * 1e9;
+            let (s, _) = real(w, LARGE).unwrap();
+            assert!(((read(&s) - w) / w).abs() < 1e-13 || w == 0.0, "{w} -> {s}");
+        }
+        for v in [0.1 + 0.2, 1.0 / 3.0, std::f64::consts::PI * 1e-7, -2.0f64.sqrt() * 1e5] {
+            let (s, _) = real(v, LARGE).unwrap();
+            assert!(((read(&s) - v) / v).abs() < 1e-13, "{v} -> {s}");
+        }
+    }
+
+    #[test]
+    fn cards_wrap_and_trim() {
+        let mut out = String::new();
+        Card::new("CHEXA").ints(1..=22u32).write(&mut out).unwrap();
+        Card::new("GRID").int(7).blank().reals([1.0, -2.5, 1.0 / 3.0]).write(&mut out).unwrap();
+        Card::new("PSHELL").int(1).int(2).real(0.5).blank().blank().write(&mut out).unwrap();
+        let want = "\
+CHEXA          1       2       3       4       5       6       7       8
+               9      10      11      12      13      14      15      16
+              17      18      19      20      21      22
+GRID*                  7                              1.            -2.5
+*       .333333333333333
+PSHELL*                1               2             0.5
+";
+        assert_eq!(out, want);
+        assert!(Card::new("X").int(123_456_789).write(&mut out).is_err(), "nine digits do not fit");
+        assert!(Card::new("X").real(f64::NAN).write(&mut out).is_err());
+    }
+}

crates/decks/src/nastran/mod.rs

@@ -0,0 +1,878 @@
+//! Nastran bulk data: the model as an input file for a linear static (SOL 101) or normal modes
+//! (SOL 103) run, written from the public Quick Reference Guide in our own words.
+//!
+//! A deck has three parts: executive control (the solution sequence), case control (one subcase
+//! per step, selecting its constraint set, load set or eigenvalue method and its output) and bulk
+//! data (the entries). What the model holds becomes:
+//!
+//! - nodes and reference points: `GRID`;
+//! - solids: `CTETRA` (tet4, tet10), `CPENTA` (wedge6), `CHEXA` (hex8, hex20) with `PSOLID`;
+//!   shells: `CTRIA3`, `CQUAD4`, `CTRIA6`, `CQUAD8` with `PSHELL`; beams: `CBEAM` with `PBEAML`;
+//!   rods: `CROD` with `PROD`. Nastran numbers a 20-node hexahedron's mid-side nodes bottom
+//!   edges, vertical edges, top edges, where the model has bottom, top, vertical
+//!   ([`grids`] reorders them); every other kind keeps the model's order;
+//! - the element formulation: `PSOLID`'s integration network and scheme (see [`solid_rule`]);
+//!   formulations Nastran's linear elements have no counterpart for are reported, not
+//!   approximated;
+//! - isotropic materials: `MAT1`; orthotropic engineering constants: `MAT9` (the stiffness matrix
+//!   in the basic axes) for solids and `MAT8` for shells lying in planes of constant Z, with the
+//!   material X axis tied to the basic X axis through a `CORD2R`;
+//! - boundary conditions: per step, one constraint set of `SPC1` (held degrees of freedom) and
+//!   `SPC` (prescribed values); in a modes step every constrained degree of freedom is held;
+//! - forces, moments and tractions: `FORCE` and `MOMENT` at the nodal values the model resolves,
+//!   so every writer and the solver share one distribution; pressure: `PLOAD4` per element face;
+//!   gravity: `GRAV`; a step with several loads combines them with `LOAD`;
+//! - kinematic coupling: `RBE2` (the reference point independent); distributing coupling: `RBE3`
+//!   (equal weights on the set's translations; the reference point dependent);
+//! - frequency steps: `EIGRL` (Lanczos) with the step's mode count and range.
+//!
+//! Field widths, continuations and number forms are in [`card`]. Model content Nastran cannot
+//! express is reported as [`Error::Unsupported`], every item at once, and nothing is written.
+//! The same model always gives the same bytes.
+
+mod card;
+
+use std::collections::{BTreeMap, BTreeSet};
+use std::fmt::Write as _;
+
+use modelcraft_model::{
+    BeamProfile, Coupling, CouplingKind, Dofs, Elastic, ElementKind, Formulation, Integration, Load, MAX_LABEL, Model, OutputVariable, Procedure,
+    Resolved, SectionKind, Step,
+};
+
+use crate::{Error, Result};
+use card::Card;
+
+/// The format's name in errors.
+pub const FORMAT: &str = "Nastran";
+
+/// Most `GRAV` entries one run may select.
+const MAX_GRAV: usize = 9;
+
+/// Grids per `SPC1` entry; longer lists are split over entries of the same set.
+const SPC1_GRIDS: usize = 1000;
+
+/// The model as a Nastran deck, with every step (a deck of some steps:
+/// [`crate::Format::write_steps`]).
+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 mut deck = Deck::new(model, r, model.steps().iter().collect());
+    deck.plan()?;
+    deck.emit()
+}
+
+/// The Nastran grid order of an element's nodes (model order in, Nastran order out).
+pub fn grids(kind: ElementKind, nodes: &[u32]) -> impl Iterator<Item = u32> + '_ {
+    order(kind).iter().filter_map(|&i| nodes.get(i).copied())
+}
+
+/// For each Nastran grid position, the position in the model's node list.
+pub fn order(kind: ElementKind) -> &'static [usize] {
+    const IDENTITY: [usize; 20] = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19];
+    // Model: 9–12 bottom edges, 13–16 top edges, 17–20 vertical edges. Nastran: bottom, vertical, top.
+    const HEX20: [usize; 20] = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15];
+    match kind {
+        ElementKind::Hex20 => &HEX20,
+        k => IDENTITY.get(..k.node_count()).unwrap_or(&IDENTITY),
+    }
+}
+
+/// The connection entry's name for an element kind.
+pub fn element_card(kind: ElementKind) -> &'static str {
+    match kind {
+        ElementKind::Rod2 => "CROD",
+        ElementKind::Beam2 => "CBEAM",
+        ElementKind::Tri3 => "CTRIA3",
+        ElementKind::Tri6 => "CTRIA6",
+        ElementKind::Quad4 => "CQUAD4",
+        ElementKind::Quad8 => "CQUAD8",
+        ElementKind::Tet4 | ElementKind::Tet10 => "CTETRA",
+        ElementKind::Wedge6 => "CPENTA",
+        ElementKind::Hex8 | ElementKind::Hex20 => "CHEXA",
+    }
+}
+
+/// `PSOLID`'s integration network (`IN`) and scheme (`ISOP`) for an element kind and
+/// formulation, or why Nastran has none.
+///
+/// - Full integration is the standard isoparametric rule: two points per direction for linear
+///   kinds, three for quadratic ones (`TWO`/`THREE`, `FULL`).
+/// - Reduced integration of a 20-node hexahedron is the 2 × 2 × 2 standard rule.
+/// - Incompatible modes of an 8-node hexahedron are Nastran's bubble functions (`BUBBLE`, with
+///   its default reduced-shear scheme).
+/// - One-point reduced integration of linear solids (with hourglass control) and the hybrid
+///   (displacement–pressure) formulation have no linear Nastran counterpart.
+pub fn solid_rule(kind: ElementKind, f: Formulation) -> std::result::Result<(Option<&'static str>, Option<&'static str>), String> {
+    if f.hybrid {
+        return Err("the hybrid (displacement–pressure) formulation".into());
+    }
+    let quadratic = kind.is_quadratic();
+    match (kind, f.incompatible_modes, f.integration) {
+        (ElementKind::Hex8, true, _) => Ok((Some("BUBBLE"), None)),
+        (_, true, _) => Err(format!("incompatible modes on {} elements", kind.name())),
+        (ElementKind::Hex20, false, Integration::Reduced) => Ok((Some("TWO"), Some("FULL"))),
+        (_, false, Integration::Reduced) => Err(format!("reduced integration of {} elements", kind.name())),
+        (_, false, Integration::Full) if quadratic => Ok((Some("THREE"), Some("FULL"))),
+        (_, false, Integration::Full) => Ok((Some("TWO"), Some("FULL"))),
+    }
+}
+
+/// The stiffness of orthotropic engineering constants for Nastran's `MAT9`, as
+/// `[[C11, C12, C13], [C22, C23, C33], [G12, G23, G13]]` (the normal block's upper triangle, then
+/// the shear moduli for xy, yz and zx in `MAT9`'s order), or `None` when it is not finite.
+pub fn orthotropic_stiffness(young: [f64; 3], poisson: [f64; 3], shear: [f64; 3]) -> Option<[[f64; 3]; 3]> {
+    let [e1, e2, e3] = young;
+    let [n12, n13, n23] = poisson;
+    // Compliance of the normal block (symmetric: ν_ij / E_i = ν_ji / E_j).
+    let (s11, s22, s33) = (1.0 / e1, 1.0 / e2, 1.0 / e3);
+    let (s12, s13, s23) = (-n12 / e1, -n13 / e1, -n23 / e2);
+    let det = s11 * (s22 * s33 - s23 * s23) - s12 * (s12 * s33 - s23 * s13) + s13 * (s12 * s23 - s22 * s13);
+    let c = [
+        [(s22 * s33 - s23 * s23) / det, (s13 * s23 - s12 * s33) / det, (s12 * s23 - s13 * s22) / det],
+        [0.0, (s11 * s33 - s13 * s13) / det, (s12 * s13 - s11 * s23) / det],
+        [0.0, 0.0, (s11 * s22 - s12 * s12) / det],
+    ];
+    let [g12, g13, g23] = shear;
+    let out = [[c[0][0], c[0][1], c[0][2]], [c[1][1], c[1][2], c[2][2]], [g12, g23, g13]];
+    out.iter().flatten().all(|v| v.is_finite()).then_some(out)
+}
+
+/// Nastran's component code of a set of degrees of freedom, such as 123 or 156.
+fn component(dofs: Dofs) -> i64 {
+    dofs.list().into_iter().fold(0, |c, d| c * 10 + i64::from(d))
+}
+
+/// The next identification number of a counter.
+fn next_id(counter: &mut u32) -> Result<u32> {
+    let id = *counter;
+    if id > MAX_LABEL {
+        return Err(Error::Limit(format!("more than {MAX_LABEL} identification numbers of one kind")));
+    }
+    *counter += 1;
+    Ok(id)
+}
+
+/// The id of `value` in `list`, adding it under a new id from `counter` when it is not there.
+fn intern<T: PartialEq>(list: &mut Vec<(u32, T)>, counter: &mut u32, value: T) -> Result<u32> {
+    if let Some((id, _)) = list.iter().find(|(_, v)| *v == value) {
+        return Ok(*id);
+    }
+    let id = next_id(counter)?;
+    list.push((id, value));
+    Ok(id)
+}
+
+/// A material's entry identification numbers.
+#[derive(Debug, Clone, Copy, Default)]
+struct MaterialIds {
+    /// `MAT1`, for isotropic materials (every element kind).
+    isotropic: Option<u32>,
+    /// `MAT9`, for orthotropic materials of solids.
+    solid: Option<u32>,
+    /// `MAT8`, for orthotropic materials of shells.
+    shell: Option<u32>,
+}
+
+/// Which element families use a material.
+#[derive(Debug, Clone, Copy, Default)]
+struct MaterialUse {
+    solid: bool,
+    shell: bool,
+    line: bool,
+}
+
+/// One property entry, and what its elements' connection entries need besides their nodes.
+enum Property {
+    Solid { mid: u32, rule: (Option<&'static str>, Option<&'static str>) },
+    Shell { mid: u32, thickness: f64, offset: Option<f64>, mcid: Option<u32> },
+    Beam { mid: u32, profile: BeamProfile, orientation: [f64; 3] },
+    Rod { mid: u32, area: f64 },
+}
+
+/// A constraint set: per node, the prescribed value of each degree of freedom.
+type Constraints = BTreeMap<u32, [Option<f64>; 6]>;
+
+/// One step's subcase: the sets it selects.
+struct Subcase<'m> {
+    step: &'m Step,
+    spc: Option<u32>,
+    load: Option<u32>,
+    method: Option<u32>,
+}
+
+struct Deck<'m> {
+    model: &'m Model,
+    r: &'m Resolved,
+    sol: u32,
+    subcases: Vec<Subcase<'m>>,
+    /// Everything the format cannot express.
+    problems: BTreeSet<String>,
+    next_sid: u32,
+    /// Material name → ids.
+    materials: BTreeMap<&'m str, MaterialIds>,
+    /// The coordinate system tying shell material axes to the basic axes, when one is needed.
+    material_csys: Option<u32>,
+    /// Properties in id order (id = position + 1), with their section's name.
+    properties: Vec<(&'m str, Property)>,
+    /// Element → property id.
+    element_pid: BTreeMap<u32, u32>,
+    /// Rigid elements: (id, coupling name, coupling).
+    rigid: Vec<(u32, &'m str, &'m Coupling)>,
+    /// Node → (dependent degrees of freedom, coupling that makes them dependent).
+    dependent: BTreeMap<u32, (Dofs, &'m str)>,
+    /// Constraint sets by id.
+    spc_sets: Vec<(u32, Constraints)>,
+    /// Load name → load set id.
+    load_sids: BTreeMap<&'m str, u32>,
+    /// Combinations: (id, load set ids).
+    combinations: Vec<(u32, Vec<u32>)>,
+}
+
+/// The solution sequence of a list of steps, or why there is none.
+fn solution(steps: &[&Step]) -> std::result::Result<u32, String> {
+    let statics = steps.iter().find(|s| matches!(s.procedure, Procedure::Static {}));
+    let modes = steps.iter().find(|s| matches!(s.procedure, Procedure::Frequency { .. }));
+    match (statics, modes) {
+        (Some(_), None) => Ok(101),
+        (None, Some(_)) => Ok(103),
+        (Some(s), Some(m)) => Err(format!(
+            "one run solves either static steps (SOL 101) or frequency steps (SOL 103), and step `{}` is static while step `{}` extracts frequencies; write them as separate decks (choose the steps to write)",
+            s.name, m.name
+        )),
+        (None, None) => Err("there is no step to write".into()),
+    }
+}
+
+impl<'m> Deck<'m> {
+    fn new(model: &'m Model, r: &'m Resolved, steps: Vec<&'m Step>) -> Deck<'m> {
+        let mut problems = BTreeSet::new();
+        let sol = solution(&steps).unwrap_or_else(|p| {
+            problems.insert(p);
+            0
+        });
+        Deck {
+            model,
+            r,
+            sol,
+            subcases: steps.into_iter().map(|step| Subcase { step, spc: None, load: None, method: None }).collect(),
+            problems,
+            next_sid: 1,
+            materials: BTreeMap::new(),
+            material_csys: None,
+            properties: Vec::new(),
+            element_pid: BTreeMap::new(),
+            rigid: Vec::new(),
+            dependent: BTreeMap::new(),
+            spc_sets: Vec::new(),
+            load_sids: BTreeMap::new(),
+            combinations: Vec::new(),
+        }
+    }
+
+    fn plan(&mut self) -> Result<()> {
+        self.plan_materials();
+        self.plan_sections()?;
+        self.plan_couplings()?;
+        self.plan_constraints()?;
+        self.plan_loads()?;
+        for sub in &mut self.subcases {
+            if matches!(sub.step.procedure, Procedure::Frequency { .. }) {
+                sub.method = Some(next_id(&mut self.next_sid)?);
+            }
+        }
+        crate::unsupported(FORMAT, std::mem::take(&mut self.problems))
+    }
+
+    /// Material ids, for the forms the sections need.
+    fn plan_materials(&mut self) {
+        let model = self.model;
+        let mut uses: BTreeMap<&str, MaterialUse> = BTreeMap::new();
+        for s in model.sections().values() {
+            let u = uses.entry(s.material.as_str()).or_default();
+            match s.kind {
+                SectionKind::Solid {} => u.solid = true,
+                SectionKind::Shell { .. } => u.shell = true,
+                SectionKind::Beam { .. } | SectionKind::Rod { .. } => u.line = true,
+            }
+        }
+        let mut next_mid = 1u32;
+        let mut take = |wanted: bool| {
+            wanted.then(|| {
+                next_mid += 1;
+                next_mid - 1
+            })
+        };
+        for (name, u) in uses {
+            let ids = match model.material(name).and_then(|m| m.elastic) {
+                Some(Elastic::Isotropic { .. }) => MaterialIds { isotropic: take(true), ..MaterialIds::default() },
+                Some(Elastic::EngineeringConstants { .. }) => {
+                    if u.line {
+                        self.problems
+                            .insert(format!("material `{name}`: orthotropic constants on beams or rods (their properties take isotropic materials)"));
+                    }
+                    MaterialIds { solid: take(u.solid), shell: take(u.shell), ..MaterialIds::default() }
+                }
+                None => MaterialIds::default(),
+            };
+            self.materials.insert(name, ids);
+        }
+    }
+
+    /// A new property for `elements`.
+    fn add_property(&mut self, section: &'m str, property: Property, elements: impl IntoIterator<Item = u32>) -> Result<()> {
+        self.properties.push((section, property));
+        let pid = u32::try_from(self.properties.len()).ok().filter(|p| *p <= MAX_LABEL).ok_or_else(|| Error::Limit("too many properties".into()))?;
+        for e in elements {
+            self.element_pid.insert(e, pid);
+        }
+        Ok(())
+    }
+
+    /// Properties and each element's property.
+    fn plan_sections(&mut self) -> Result<()> {
+        let (model, r) = (self.model, self.r);
+        for (sname, sec) in model.sections() {
+            let ids = self.materials.get(sec.material.as_str()).copied().unwrap_or_default();
+            let elements = r.sets.get(&sec.set).map(|s| s.elements.as_slice()).unwrap_or_default();
+            let kind_of = |e: u32| model.element(e).map(|el| el.kind);
+            match sec.kind {
+                SectionKind::Solid {} => {
+                    let Some(mid) = ids.isotropic.or(ids.solid) else { continue };
+                    // One property per integration rule the section's element kinds need.
+                    let kinds: BTreeSet<ElementKind> = elements.iter().filter_map(|e| kind_of(*e)).collect();
+                    let mut rules = BTreeMap::new();
+                    for kind in kinds {
+                        match solid_rule(kind, sec.formulation) {
+                            Ok(rule) => {
+                                rules.entry(rule).or_insert_with(Vec::new).push(kind);
+                            }
+                            Err(why) => {
+                                self.problems.insert(format!("section `{sname}`: {why} (Nastran's linear solids have no such formulation)"));
+                            }
+                        }
+                    }
+                    for (rule, kinds) in rules {
+                        let members = elements.iter().copied().filter(|e| kind_of(*e).is_some_and(|k| kinds.contains(&k)));
+                        self.add_property(sname, Property::Solid { mid, rule }, members)?;
+                    }
+                }
+                SectionKind::Shell { thickness, offset, points } => {
+                    if sec.formulation.integration == Integration::Reduced {
+                        self.problems.insert(format!("section `{sname}`: reduced integration of shells (Nastran's shells have one formulation)"));
+                    }
+                    if points < 2 {
+                        self.problems.insert(format!(
+                            "section `{sname}`: one integration point through the thickness (no bending stiffness; Nastran integrates a linear section exactly)"
+                        ));
+                    }
+                    let (mid, mcid) = match ids {
+                        MaterialIds { isotropic: Some(m), .. } => (m, None),
+                        MaterialIds { shell: Some(m), .. } => {
+                            if !elements.iter().all(|e| model.element(*e).is_some_and(|el| self.in_constant_z(&el.nodes, el.kind))) {
+                                self.problems.insert(format!(
+                                    "section `{sname}`: orthotropic shells that do not lie in planes of constant Z (the material axes are the global ones)"
+                                ));
+                            }
+                            (m, Some(*self.material_csys.get_or_insert(1)))
+                        }
+                        _ => continue,
+                    };
+                    // The nodes lie on the reference surface; Nastran's offset runs from them to
+                    // the mid-surface along the element normal.
+                    let offset = (offset != 0.0).then_some(-offset * thickness);
+                    self.add_property(sname, Property::Shell { mid, thickness, offset, mcid }, elements.iter().copied())?;
+                }
+                SectionKind::Beam { profile, orientation } => {
+                    let Some(mid) = ids.isotropic else { continue };
+                    self.add_property(sname, Property::Beam { mid, profile, orientation }, elements.iter().copied())?;
+                }
+                SectionKind::Rod { area } => {
+                    let Some(mid) = ids.isotropic else { continue };
+                    self.add_property(sname, Property::Rod { mid, area }, elements.iter().copied())?;
+                }
+            }
+        }
+        Ok(())
+    }
+
+    /// Whether a shell element's corners lie in a plane of constant Z.
+    fn in_constant_z(&self, nodes: &[u32], kind: ElementKind) -> bool {
+        let corners = nodes.get(..kind.corner_count()).unwrap_or_default();
+        let mut z = corners.iter().map(|n| self.model.node(*n).map(|n| n.xyz[2]));
+        let first = z.next().flatten();
+        first.is_some() && z.all(|v| v == first)
+    }
+
+    /// Rigid element ids and the degrees of freedom couplings make dependent.
+    fn plan_couplings(&mut self) -> Result<()> {
+        let (model, r) = (self.model, self.r);
+        let mut next = model.elements().keys().next_back().copied().unwrap_or(0);
+        let mut overlaps: BTreeMap<(&str, &str), BTreeSet<u32>> = BTreeMap::new();
+        for (name, c) in model.couplings() {
+            next = next
+                .checked_add(1)
+                .filter(|n| *n <= MAX_LABEL)
+                .ok_or_else(|| Error::Limit(format!("coupling `{name}`: no element identification number left below {}", MAX_LABEL as u64 + 1)))?;
+            self.rigid.push((next, name.as_str(), c));
+            let dependents: &[u32] = match c.kind {
+                CouplingKind::Kinematic => r.sets.get(&c.set).map(|s| s.nodes.as_slice()).unwrap_or_default(),
+                CouplingKind::Distributing => r.reference_point_nodes.get(&c.reference_point).map(std::slice::from_ref).unwrap_or_default(),
+            };
+            for &node in dependents {
+                match self.dependent.get(&node).copied() {
+                    Some((had, other)) if had.bits() & c.dofs.bits() != 0 => {
+                        overlaps.entry((other, name.as_str())).or_default().insert(node);
+                    }
+                    Some((had, other)) => {
+                        self.dependent.insert(node, (Dofs::from_bits(had.bits() | c.dofs.bits()), other));
+                    }
+                    None => {
+                        self.dependent.insert(node, (c.dofs, name));
+                    }
+                }
+            }
+        }
+        for ((a, b), nodes) in overlaps {
+            self.problems.insert(format!(
+                "couplings `{a}` and `{b}` both make nodes {} follow another node (Nastran makes a degree of freedom dependent once)",
+                labels(&nodes)
+            ));
+        }
+        Ok(())
+    }
+
+    /// One constraint set per distinct step constraint, checked against coupled nodes.
+    fn plan_constraints(&mut self) -> Result<()> {
+        let (model, r) = (self.model, self.r);
+        // Conflicts, one problem per pair of definitions: the nodes they meet on.
+        let mut clashes: BTreeMap<(&str, &str, &str), BTreeSet<u32>> = BTreeMap::new();
+        let mut on_dependent: BTreeMap<(&str, &str), BTreeSet<u32>> = BTreeMap::new();
+        for sub in &mut self.subcases {
+            let step = sub.step;
+            let modes = matches!(step.procedure, Procedure::Frequency { .. });
+            let mut set: Constraints = BTreeMap::new();
+            let mut by: BTreeMap<(u32, u8), &str> = BTreeMap::new();
+            for name in &step.boundary_conditions {
+                let Some(bc) = model.boundary_condition(name) else { continue };
+                let dofs = bc.dofs();
+                // In a modes step every constrained degree of freedom is held.
+                let values: Vec<(u8, f64)> = dofs.list().into_iter().filter_map(|d| Some((d, if modes { 0.0 } else { bc.value(d)? }))).collect();
+                for &node in r.target_nodes(bc.target()) {
+                    if let Some((dep, coupling)) = self.dependent.get(&node)
+                        && dep.bits() & dofs.bits() != 0
+                    {
+                        on_dependent.entry((name.as_str(), coupling)).or_default().insert(node);
+                    }
+                    let slots = set.entry(node).or_default();
+                    for &(dof, value) in &values {
+                        let Some(slot) = usize::from(dof).checked_sub(1).and_then(|i| slots.get_mut(i)) else { continue };
+                        match *slot {
+                            Some(v) if v != value => {
+                                let other = by.get(&(node, dof)).copied().unwrap_or_default();
+                                clashes.entry((step.name.as_str(), other, name.as_str())).or_default().insert(node);
+                            }
+                            _ => {
+                                *slot = Some(value);
+                                by.insert((node, dof), name);
+                            }
+                        }
+                    }
+                }
+            }
+            if !set.is_empty() {
+                sub.spc = Some(intern(&mut self.spc_sets, &mut self.next_sid, set)?);
+            }
+        }
+        for ((step, a, b), nodes) in clashes {
+            self.problems.insert(format!(
+                "step `{step}`: boundary conditions `{a}` and `{b}` prescribe different values at nodes {} (Nastran takes one value per degree of freedom)",
+                labels(&nodes)
+            ));
+        }
+        for ((bc, coupling), nodes) in on_dependent {
+            self.problems.insert(format!(
+                "boundary condition `{bc}` constrains nodes {} whose motion coupling `{coupling}` prescribes (Nastran cannot constrain a dependent degree of freedom)",
+                labels(&nodes)
+            ));
+        }
+        Ok(())
+    }
+
+    /// Load set ids for the static steps' loads, and their combinations.
+    fn plan_loads(&mut self) -> Result<()> {
+        let (model, r) = (self.model, self.r);
+        let mut gravity = 0usize;
+        for sub in &mut self.subcases {
+            if !matches!(sub.step.procedure, Procedure::Static {}) {
+                continue;
+            }
+            let mut sids = Vec::new();
+            for name in &sub.step.loads {
+                if let Some(sid) = self.load_sids.get(name.as_str()) {
+                    sids.push(*sid);
+                    continue;
+                }
+                let Some(load) = model.load(name) else { continue };
+                match load {
+                    Load::BodyForce { .. } => {
+                        self.problems.insert(format!("load `{name}`: a body force (Nastran has no force-per-volume load for these elements)"));
+                    }
+                    Load::Gravity { set: Some(s), .. } if r.sets.get(s).map(|s| s.elements.len()) != Some(model.elements().len()) => {
+                        self.problems.insert(format!("load `{name}`: gravity on part of the model (set `{s}`); Nastran's gravity acts on all of it"));
+                    }
+                    Load::Gravity { .. } => gravity += 1,
+                    _ => {}
+                }
+                let sid = next_id(&mut self.next_sid)?;
+                self.load_sids.insert(name.as_str(), sid);
+                sids.push(sid);
+            }
+            sub.load = match sids.as_slice() {
+                [] => None,
+                [one] => Some(*one),
+                _ => Some(intern(&mut self.combinations, &mut self.next_sid, sids)?),
+            };
+        }
+        if gravity > MAX_GRAV {
+            self.problems.insert(format!("{gravity} gravity loads (a run selects at most {MAX_GRAV})"));
+        }
+        Ok(())
+    }
+
+    fn emit(&self) -> Result<String> {
+        let mut out = String::new();
+        out.push_str("$ Nastran bulk data written by ModelCraft\n");
+        let _ = writeln!(out, "$ units: {}", self.model.units().label());
+        let _ = writeln!(out, "SOL {}\nCEND\nTITLE = ModelCraft model\nECHO = NONE", self.sol);
+        for (i, sub) in self.subcases.iter().enumerate() {
+            let _ = writeln!(out, "SUBCASE {}", i + 1);
+            let title: String = format!("  SUBTITLE = {}", sub.step.name).chars().take(72).collect();
+            let _ = writeln!(out, "{title}");
+            for (key, sid) in [("SPC", sub.spc), ("LOAD", sub.load), ("METHOD", sub.method)] {
+                if let Some(sid) = sid {
+                    let _ = writeln!(out, "  {key} = {sid}");
+                }
+            }
+            for v in &sub.step.output {
+                let request = match v {
+                    OutputVariable::Displacement => "DISPLACEMENT",
+                    OutputVariable::ReactionForce => "SPCFORCES",
+                    OutputVariable::Stress => "STRESS",
+                    OutputVariable::Strain => "STRAIN",
+                };
+                let _ = writeln!(out, "  {request} = ALL");
+            }
+        }
+        out.push_str("BEGIN BULK\n");
+        Card::new("PARAM").text("POST").int(-1).write(&mut out)?;
+        Card::new("PARAM").text("AUTOSPC").text("YES").write(&mut out)?;
+        if self.sol == 103 {
+            // Consistent (coupled) mass, as the other formats use by default.
+            Card::new("PARAM").text("COUPMASS").int(1).write(&mut out)?;
+        }
+        self.emit_mesh(&mut out)?;
+        self.emit_constraints(&mut out)?;
+        self.emit_loads(&mut out)?;
+        if self.subcases.iter().any(|s| s.method.is_some()) {
+            comment(&mut out, "eigenvalue extraction");
+        }
+        for sub in &self.subcases {
+            if let (Some(sid), Procedure::Frequency { modes, min_hz, max_hz }) = (sub.method, sub.step.procedure) {
+                Card::new("EIGRL").int(sid).opt(min_hz).opt(max_hz).int(modes).write(&mut out)?;
+            }
+        }
+        out.push_str("ENDDATA\n");
+        Ok(out)
+    }
+
+    /// Grids, materials, properties, elements and couplings.
+    fn emit_mesh(&self, out: &mut String) -> Result<()> {
+        comment(out, "nodes");
+        for (label, node) in self.model.nodes() {
+            Card::new("GRID").int(*label).blank().reals(node.xyz).write(out)?;
+        }
+        if !self.model.reference_points().is_empty() {
+            comment(out, "reference points");
+        }
+        for (name, rp) in self.model.reference_points() {
+            comment(out, &format!("reference point {name}: GRID {}", rp.label));
+            Card::new("GRID").int(rp.label).blank().reals(rp.xyz).write(out)?;
+        }
+        if let Some(cid) = self.material_csys {
+            comment(out, "material axes of orthotropic shells: the basic axes");
+            Card::new("CORD2R").int(cid).blank().reals([0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0]).write(out)?;
+        }
+
+        comment(out, "materials");
+        for (name, ids) in &self.materials {
+            let Some(m) = self.model.material(name) else { continue };
+            let rho = m.density;
+            let card = match (m.elastic, ids) {
+                (Some(Elastic::Isotropic { young, poisson }), MaterialIds { isotropic: Some(mid), .. }) => {
+                    Some(Card::new("MAT1").int(*mid).real(young).blank().real(poisson).opt(rho))
+                }
+                _ => None,
+            };
+            if let (Some(Elastic::EngineeringConstants { young, poisson, shear }), Some(mid)) = (m.elastic, ids.solid) {
+                let [[c11, c12, c13], [c22, c23, c33], [g12, g23, g13]] = orthotropic_stiffness(young, poisson, shear)
+                    .ok_or_else(|| Error::Limit(format!("material `{name}`: its stiffness is not a finite number")))?;
+                let card = Card::new("MAT9")
+                    .int(mid)
+                    .reals([c11, c12, c13, 0.0, 0.0, 0.0, c22, c23, 0.0, 0.0, 0.0, c33, 0.0, 0.0, 0.0, g12, 0.0, 0.0, g23, 0.0, g13])
+                    .opt(rho);
+                write_named(out, &format!("material {name}"), mid, card)?;
+            }
+            if let (Some(Elastic::EngineeringConstants { young: [e1, e2, _], poisson: [nu12, ..], shear: [g12, g13, g23] }), Some(mid)) =
+                (m.elastic, ids.shell)
+            {
+                write_named(out, &format!("material {name}"), mid, Card::new("MAT8").int(mid).reals([e1, e2, nu12, g12, g13, g23]).opt(rho))?;
+            }
+            if let (Some(card), Some(mid)) = (card, ids.isotropic) {
+                write_named(out, &format!("material {name}"), mid, card)?;
+            }
+        }
+
+        comment(out, "properties");
+        for (pid, (section, p)) in (1u32..).zip(&self.properties) {
+            let card = match p {
+                Property::Solid { mid, rule: (network, scheme) } => Card::new("PSOLID").int(pid).int(*mid).blank().opt(*network).blank().opt(*scheme),
+                Property::Shell { mid, thickness, .. } => Card::new("PSHELL").int(pid).int(*mid).real(*thickness).int(*mid).blank().int(*mid),
+                Property::Beam { mid, profile, .. } => {
+                    let (shape, dims) = beam_dimensions(*profile);
+                    Card::new("PBEAML").int(pid).int(*mid).blank().text(shape).end_line().reals(dims)
+                }
+                Property::Rod { mid, area } => Card::new("PROD").int(pid).int(*mid).real(*area),
+            };
+            write_named(out, &format!("section {section}"), pid, card)?;
+        }
+
+        comment(out, "elements");
+        for (label, el) in self.model.elements() {
+            let Some(pid) = self.element_pid.get(label).copied() else { continue };
+            let Some((_, property)) = usize::try_from(pid).ok().and_then(|p| self.properties.get(p - 1)) else { continue };
+            let card = Card::new(element_card(el.kind)).int(*label).int(pid).ints(grids(el.kind, &el.nodes));
+            let card = match property {
+                Property::Beam { orientation, .. } => card.reals(self.beam_vector(&el.nodes, *orientation)),
+                // CQUAD8 has four corner thicknesses before the material axes and the offset.
+                Property::Shell { offset, mcid, .. } if el.kind == ElementKind::Quad8 => card.blank().blank().blank().blank().opt(*mcid).opt(*offset),
+                Property::Shell { offset, mcid, .. } => card.opt(*mcid).opt(*offset),
+                Property::Solid { .. } | Property::Rod { .. } => card,
+            };
+            card.write(out)?;
+        }
+
+        if !self.rigid.is_empty() {
+            comment(out, "couplings");
+        }
+        for (eid, name, c) in &self.rigid {
+            let rp = self.r.reference_point_nodes.get(&c.reference_point).copied().unwrap_or_default();
+            let nodes = self.r.sets.get(&c.set).map(|s| s.nodes.as_slice()).unwrap_or_default();
+            let card = match c.kind {
+                CouplingKind::Kinematic => Card::new("RBE2").int(*eid).int(rp).int(component(c.dofs)).ints(nodes.iter().copied()),
+                // Equal weights on the set's translations.
+                CouplingKind::Distributing => {
+                    Card::new("RBE3").int(*eid).blank().int(rp).int(component(c.dofs)).real(1.0).int(123).ints(nodes.iter().copied())
+                }
+            };
+            write_named(out, &format!("coupling {name}"), *eid, card)?;
+        }
+        Ok(())
+    }
+
+    /// Constraint sets: held degrees of freedom by `SPC1`, prescribed values by `SPC`.
+    fn emit_constraints(&self, out: &mut String) -> Result<()> {
+        if !self.spc_sets.is_empty() {
+            comment(out, "boundary conditions");
+        }
+        for (sid, set) in &self.spc_sets {
+            for sub in self.subcases.iter().filter(|s| s.spc == Some(*sid)) {
+                comment(out, &format!("step {}: SPC {sid}", sub.step.name));
+            }
+            let mut held: BTreeMap<u8, Vec<u32>> = BTreeMap::new();
+            let mut prescribed = Vec::new();
+            for (node, values) in set {
+                let mut zero = 0u8;
+                // Nonzero values: one entry per value, with every degree of freedom it applies to.
+                let mut by_value: Vec<(f64, u8)> = Vec::new();
+                for (bit, v) in values.iter().enumerate().filter_map(|(i, v)| Some((1u8 << i, (*v)?))) {
+                    match by_value.iter_mut().find(|(w, _)| *w == v) {
+                        _ if v == 0.0 => zero |= bit,
+                        Some((_, m)) => *m |= bit,
+                        None => by_value.push((v, bit)),
+                    }
+                }
+                if zero != 0 {
+                    held.entry(zero).or_default().push(*node);
+                }
+                prescribed.extend(by_value.into_iter().map(|(v, m)| (*node, m, v)));
+            }
+            for (mask, nodes) in &held {
+                for chunk in nodes.chunks(SPC1_GRIDS) {
+                    Card::new("SPC1").int(*sid).int(component(Dofs::from_bits(*mask))).ints(chunk.iter().copied()).write(out)?;
+                }
+            }
+            for (node, mask, value) in prescribed {
+                Card::new("SPC").int(*sid).int(node).int(component(Dofs::from_bits(mask))).real(value).write(out)?;
+            }
+        }
+        Ok(())
+    }
+
+    /// Load sets and their combinations.
+    fn emit_loads(&self, out: &mut String) -> Result<()> {
+        if !self.load_sids.is_empty() {
+            comment(out, "loads");
+        }
+        for (name, sid) in &self.load_sids {
+            let Some(load) = self.model.load(name) else { continue };
+            let entry = |card: Card, out: &mut String, first: &mut bool| -> Result<()> {
+                if std::mem::take(first) {
+                    comment(out, &format!("load {name}: {} {sid}", card.name()));
+                }
+                card.write(out)
+            };
+            let mut first = true;
+            match load {
+                Load::Force { .. } | Load::Traction { .. } | Load::Moment { .. } => {
+                    let card = if matches!(load, Load::Moment { .. }) { "MOMENT" } else { "FORCE" };
+                    let values = self.r.nodal_loads.get(*name).map(Vec::as_slice).unwrap_or_default();
+                    for (node, v) in values.iter().filter(|(_, v)| v.iter().any(|c| *c != 0.0)) {
+                        entry(Card::new(card).int(*sid).int(*node).blank().real(1.0).reals(*v), out, &mut first)?;
+                    }
+                }
+                Load::Pressure { set, magnitude } => {
+                    for f in self.r.sets.get(set).map(|s| s.faces.as_slice()).unwrap_or_default() {
+                        if let Some(card) = self.pressure(*sid, *f, *magnitude) {
+                            entry(card, out, &mut first)?;
+                        }
+                    }
+                }
+                Load::Gravity { acceleration, .. } => entry(Card::new("GRAV").int(*sid).blank().real(1.0).reals(*acceleration), out, &mut first)?,
+                Load::BodyForce { .. } => {}
+            }
+        }
+        for (sid, loads) in &self.combinations {
+            let mut card = Card::new("LOAD").int(*sid).real(1.0);
+            for l in loads {
+                card = card.real(1.0).int(*l);
+            }
+            card.write(out)?;
+        }
+        Ok(())
+    }
+
+    /// A `PLOAD4` for one element face, positive pushing into the element.
+    fn pressure(&self, sid: u32, f: modelcraft_model::ElementFace, magnitude: f64) -> Option<Card> {
+        let el = self.model.element(f.element)?;
+        let card = Card::new("PLOAD4").int(sid).int(f.element);
+        if el.kind.dimension() == 2 {
+            // Nastran pushes along the element normal; face 1 is the side the normal points to.
+            let p = if f.face == 1 { -magnitude } else { magnitude };
+            return Some(card.real(p));
+        }
+        let corners = el.kind.corner_count();
+        let local = el.kind.faces().get(usize::from(f.face).checked_sub(1)?)?;
+        let face_corners = local.get(..local.iter().take_while(|&&i| i < corners).count())?;
+        let node = |i: usize| el.nodes.get(i).copied();
+        let card = card.real(magnitude).blank().blank().blank().int(node(*face_corners.first()?)?);
+        Some(match (el.kind, face_corners.len()) {
+            // Tetrahedra: the corner off the face.
+            (ElementKind::Tet4 | ElementKind::Tet10, _) => card.int(node((0..corners).find(|i| !face_corners.contains(i))?)?),
+            // Quadrilateral faces: the corner diagonally opposite G1.
+            (_, 4) => card.int(node(*face_corners.get(2)?)?),
+            // Triangular faces of wedges: G3 stays blank.
+            _ => card,
+        })
+    }
+
+    /// The `CBEAM` orientation vector: the profile's local 1 direction (the section's
+    /// orientation) is Nastran's element z axis, so v = axis × orientation.
+    fn beam_vector(&self, nodes: &[u32], o: [f64; 3]) -> [f64; 3] {
+        let p = |i: usize| nodes.get(i).and_then(|n| self.model.node(*n)).map(|n| n.xyz).unwrap_or_default();
+        let (a, b) = (p(0), p(1));
+        let t = [b[0] - a[0], b[1] - a[1], b[2] - a[2]];
+        [t[1] * o[2] - t[2] * o[1], t[2] * o[0] - t[0] * o[2], t[0] * o[1] - t[1] * o[0]]
+    }
+}
+
+/// `PBEAML`'s shape and dimensions. Nastran's sections have their height along the element y
+/// axis and their width (flanges) along z, which [`Deck::beam_vector`] makes the profile's
+/// local 1 direction.
+fn beam_dimensions(p: BeamProfile) -> (&'static str, Vec<f64>) {
+    match p {
+        BeamProfile::Rectangle { width, height } => ("BAR", vec![width, height]),
+        BeamProfile::Circle { radius } => ("ROD", vec![radius]),
+        BeamProfile::Pipe { radius, thickness } if radius - thickness > 0.0 => ("TUBE", vec![radius, radius - thickness]),
+        BeamProfile::Pipe { radius, .. } => ("ROD", vec![radius]),
+        BeamProfile::I { height, width, web, flange } => ("I", vec![height, width, width, web, flange, flange]),
+        BeamProfile::Box { width, height, thickness } => ("BOX", vec![width, height, thickness, thickness]),
+    }
+}
+
+/// An entry preceded by a comment naming what it holds (`material Steel: MAT1 1`).
+fn write_named(out: &mut String, what: &str, id: u32, card: Card) -> Result<()> {
+    comment(out, &format!("{what}: {} {id}", card.name()));
+    card.write(out)
+}
+
+/// The first few labels of a list, and how many more there are.
+fn labels(nodes: &BTreeSet<u32>) -> String {
+    const SHOWN: usize = 5;
+    let shown: Vec<String> = nodes.iter().take(SHOWN).map(u32::to_string).collect();
+    match nodes.len().saturating_sub(SHOWN) {
+        0 => shown.join(", "),
+        more => format!("{}, and {more} more", shown.join(", ")),
+    }
+}
+
+/// A comment line, at most 80 columns.
+fn comment(out: &mut String, text: &str) {
+    out.extend(format!("$ {text}").chars().take(80));
+    out.push('\n');
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn orthotropic_stiffness_of_isotropic_constants_is_lames() {
+        let (e, nu) = (210_000.0, 0.3);
+        let g = e / (2.0 * (1.0 + nu));
+        let c = orthotropic_stiffness([e; 3], [nu; 3], [g; 3]).unwrap();
+        let k = e / ((1.0 + nu) * (1.0 - 2.0 * nu));
+        let close = |a: f64, b: f64| (a - b).abs() < 1e-9 * b.abs();
+        assert!(close(c[0][0], k * (1.0 - nu)) && close(c[1][0], k * (1.0 - nu)) && close(c[1][2], k * (1.0 - nu)));
+        assert!(close(c[0][1], k * nu) && close(c[0][2], k * nu) && close(c[1][1], k * nu));
+        assert_eq!(c[2], [g; 3]);
+    }
+
+    #[test]
+    fn hex20_mid_side_nodes_swap_blocks() {
+        let nodes: Vec<u32> = (1..=20).collect();
+        let g: Vec<u32> = grids(ElementKind::Hex20, &nodes).collect();
+        assert_eq!(&g[..12], &nodes[..12]);
+        assert_eq!(&g[12..16], &[17, 18, 19, 20]);
+        assert_eq!(&g[16..], &[13, 14, 15, 16]);
+        for k in ElementKind::ALL {
+            let mut sorted = order(k).to_vec();
+            sorted.sort_unstable();
+            assert_eq!(sorted, (0..k.node_count()).collect::<Vec<_>>(), "{k:?} is a permutation");
+        }
+    }
+
+    #[test]
+    fn formulations_map_or_are_named() {
+        let f = |integration, incompatible_modes, hybrid| Formulation { integration, incompatible_modes, hybrid };
+        assert_eq!(solid_rule(ElementKind::Hex8, f(Integration::Full, false, false)), Ok((Some("TWO"), Some("FULL"))));
+        assert_eq!(solid_rule(ElementKind::Hex8, f(Integration::Full, true, false)), Ok((Some("BUBBLE"), None)));
+        assert_eq!(solid_rule(ElementKind::Hex20, f(Integration::Reduced, false, false)), Ok((Some("TWO"), Some("FULL"))));
+        assert_eq!(solid_rule(ElementKind::Tet10, f(Integration::Full, false, false)), Ok((Some("THREE"), Some("FULL"))));
+        assert!(solid_rule(ElementKind::Hex8, f(Integration::Reduced, false, false)).is_err());
+        assert!(solid_rule(ElementKind::Tet4, f(Integration::Full, false, true)).is_err());
+    }
+}

crates/decks/tests/golden/nastran/cantilever.bdf

@@ -0,0 +1,78 @@
+$ Nastran bulk data written by ModelCraft
+$ units: mm-t-s-N-MPa
+SOL 101
+CEND
+TITLE = ModelCraft model
+ECHO = NONE
+SUBCASE 1
+  SUBTITLE = Bend
+  SPC = 1
+  LOAD = 2
+  DISPLACEMENT = ALL
+  SPCFORCES = ALL
+  STRESS = ALL
+  STRAIN = ALL
+BEGIN BULK
+PARAM       POST      -1
+PARAM    AUTOSPC     YES
+$ nodes
+GRID*                  1                              0.              0.
+*                     0.
+GRID*                  2                              0.              0.
+*                     1.
+GRID*                  3                              0.              1.
+*                     0.
+GRID*                  4                              0.              1.
+*                     1.
+GRID*                  5                              1.              0.
+*                     0.
+GRID*                  6                              1.              0.
+*                     1.
+GRID*                  7                              1.              1.
+*                     0.
+GRID*                  8                              1.              1.
+*                     1.
+GRID*                  9                              2.              0.
+*                     0.
+GRID*                 10                              2.              0.
+*                     1.
+GRID*                 11                              2.              1.
+*                     0.
+GRID*                 12                              2.              1.
+*                     1.
+GRID*                 13                              3.              0.
+*                     0.
+GRID*                 14                              3.              0.
+*                     1.
+GRID*                 15                              3.              1.
+*                     0.
+GRID*                 16                              3.              1.
+*                     1.
+$ materials
+$ material Steel: MAT1 1
+MAT1*                  1         210000.                             0.3
+*                7.85E-9
+$ properties
+$ section bar: PSOLID 1
+PSOLID         1       1             TWO            FULL
+$ elements
+CHEXA          1       1       1       5       7       3       2       6
+               8       4
+CHEXA          2       1       5       9      11       7       6      10
+              12       8
+CHEXA          3       1       9      13      15      11      10      14
+              16      12
+$ boundary conditions
+$ step Bend: SPC 1
+SPC1           1  123456       1       2       3       4
+$ loads
+$ load tip: FORCE 2
+FORCE*                 2              13                              1.
+*                     0.              0.            -25.
+FORCE*                 2              14                              1.
+*                     0.              0.            -25.
+FORCE*                 2              15                              1.
+*                     0.              0.            -25.
+FORCE*                 2              16                              1.
+*                     0.              0.            -25.
+ENDDATA

crates/decks/tests/golden/nastran/modes.bdf

@@ -0,0 +1,82 @@
+$ Nastran bulk data written by ModelCraft
+$ units: unspecified
+SOL 103
+CEND
+TITLE = ModelCraft model
+ECHO = NONE
+SUBCASE 1
+  SUBTITLE = Modes
+  SPC = 1
+  METHOD = 2
+  DISPLACEMENT = ALL
+BEGIN BULK
+PARAM       POST      -1
+PARAM    AUTOSPC     YES
+PARAM   COUPMASS       1
+$ nodes
+GRID*                 81                           -0.05            1.45
+*                   0.45
+GRID*                 82                            2.15             1.5
+*                   0.65
+GRID*                 83                              2.            0.15
+*                   0.55
+GRID*                 84                            0.05              0.
+*                    0.5
+GRID*                 85                            0.05             0.7
+*                   0.95
+GRID*                 86                            1.05             1.5
+*                   1.05
+GRID*                 87                            2.05             0.9
+*                   1.15
+GRID*                 88                            1.05             0.1
+*                   1.05
+GRID*                 89                           -0.05            0.75
+*                     0.
+GRID*                 90                            1.05            1.45
+*                   0.05
+GRID*                 91                             2.1            0.75
+*                   0.05
+GRID*                 92                              1.            0.05
+*                     0.
+GRID*                 93                              0.             1.4
+*                    0.9
+GRID*                 94                             2.1             1.6
+*                    1.2
+GRID*                 95                              2.             0.2
+*                    1.1
+GRID*                 96                             0.1              0.
+*                     1.
+GRID*                 97                            -0.1             1.5
+*                     0.
+GRID*                 98                             2.2             1.4
+*                    0.1
+GRID*                 99                              2.             0.1
+*                     0.
+GRID*                100                              0.              0.
+*                     0.
+$ reference points
+$ reference point RP: GRID 99999999
+GRID*           99999999                              1.            0.75
+*                     2.
+$ materials
+$ material Steel: MAT1 1
+MAT1*                  1         210000.                             0.3
+*                7.85E-9
+$ properties
+$ section s: PSOLID 1
+PSOLID         1       1           THREE            FULL
+$ elements
+CHEXA          7       1     100      99      98      97      96      95
+              94      93      92      91      90      89      84      83
+              82      81      88      87      86      85
+$ couplings
+$ coupling lid: RBE2 8
+RBE2           899999999  123456      85      86      87      88      93
+              94      95      96
+$ boundary conditions
+$ step Modes: SPC 1
+SPC1           1  123456      89      90      91      92      97      98
+              99     100
+$ eigenvalue extraction
+EIGRL*                 2                           5000.               6
+ENDDATA

crates/decks/tests/golden/nastran/plate.bdf

@@ -0,0 +1,64 @@
+$ Nastran bulk data written by ModelCraft
+$ units: mm-t-s-N-MPa
+SOL 101
+CEND
+TITLE = ModelCraft model
+ECHO = NONE
+SUBCASE 1
+  SUBTITLE = Sag
+  SPC = 1
+  LOAD = 4
+  DISPLACEMENT = ALL
+  SPCFORCES = ALL
+  STRESS = ALL
+  STRAIN = ALL
+BEGIN BULK
+PARAM       POST      -1
+PARAM    AUTOSPC     YES
+$ nodes
+GRID*                  1                              0.              0.
+*                     0.
+GRID*                  2                              5.              0.
+*                     0.
+GRID*                  3                             10.              0.
+*                     0.
+GRID*                  4                              0.              5.
+*                     0.
+GRID*                  5                              5.              5.
+*                     0.
+GRID*                  6                             10.              5.
+*                     0.
+GRID*                  7                              0.             10.
+*                     0.
+GRID*                  8                              5.             10.
+*                     0.
+GRID*                  9                             10.             10.
+*                     0.
+$ materials
+$ material Aluminium: MAT1 1
+MAT1*                  1          70000.                            0.33
+*                 2.7E-9
+$ properties
+$ section skin: PSHELL 1
+PSHELL*                1               1              2.               1
+*                                      1
+$ elements
+CQUAD4         1       1       1       2       5       4
+CQUAD4         2       1       2       3       6       5
+CQUAD4         3       1       4       5       8       7
+CQUAD4         4       1       5       6       9       8
+$ boundary conditions
+$ step Sag: SPC 1
+SPC1           1  123456       1       4       7
+$ loads
+$ load press: PLOAD4 2
+PLOAD4*                2               1           -0.01
+PLOAD4*                2               2           -0.01
+PLOAD4*                2               3           -0.01
+PLOAD4*                2               4           -0.01
+$ load weight: GRAV 3
+GRAV*                  3                              1.              0.
+*                     0.          -9810.
+LOAD*                  4              1.              1.               2
+*                     1.               3
+ENDDATA

crates/decks/tests/nastran.rs

@@ -0,0 +1,732 @@
+//! The Nastran writer on our own small models: golden files, a minimal re-reader of our own
+//! output (grids, elements, constraint sets, case control), node order checked against Nastran's
+//! own element conventions by Jacobian sign and mid-side positions, and what is unsupported.
+//!
+//! Golden files are rewritten with `MODELCRAFT_BLESS=1 cargo test -p modelcraft-decks`.
+
+use std::collections::BTreeMap;
+
+use modelcraft_decks::nastran;
+use modelcraft_decks::{Error, with_steps};
+use modelcraft_model::{
+    BeamProfile, BoundaryCondition, Coupling, CouplingKind, Distribution, Dofs, Elastic, ElementFace, ElementKind, Formulation, Integration, Load,
+    Material, Model, Procedure, ReferencePoint, Section, SectionKind, Set, Step, Target, UnitSystem,
+};
+
+// ---------------------------------------------------------------------------------------------
+// A minimal reader of our own output.
+
+#[derive(Debug, Default)]
+struct Read {
+    sol: u32,
+    /// Per subcase: its case control lines as `KEY = value`.
+    subcases: Vec<BTreeMap<String, String>>,
+    /// Bulk entries: name and data fields (trimmed, blanks as "").
+    cards: Vec<(String, Vec<String>)>,
+}
+
+fn read(deck: &str) -> Read {
+    let mut r = Read::default();
+    let mut lines = deck.lines();
+    for line in lines.by_ref() {
+        let l = line.trim();
+        if l.starts_with('$') {
+            continue;
+        }
+        if let Some(s) = l.strip_prefix("SOL ") {
+            r.sol = s.parse().unwrap();
+        } else if l.starts_with("SUBCASE") {
+            r.subcases.push(BTreeMap::new());
+        } else if l == "BEGIN BULK" {
+            break;
+        } else if let (Some((k, v)), Some(sub)) = (l.split_once('='), r.subcases.last_mut()) {
+            sub.insert(k.trim().to_string(), v.trim().to_string());
+        }
+    }
+    for line in lines {
+        assert!(line.len() <= 80, "{line}");
+        assert!(!line.contains('\t'));
+        if line.starts_with('$') || line == "ENDDATA" {
+            continue;
+        }
+        let (first, rest) = line.split_at(8.min(line.len()));
+        let large = first.trim_end().ends_with('*');
+        let width = if large { 16 } else { 8 };
+        let per_line = if large { 4 } else { 8 };
+        // Every physical line holds its full count of fields; trailing blanks are trimmed away.
+        let mut fields: Vec<String> = rest.as_bytes().chunks(width).map(|c| String::from_utf8_lossy(c).trim().to_string()).collect();
+        fields.resize(per_line, String::new());
+        let continuation = first.trim().is_empty() || first.trim() == "*";
+        if continuation {
+            r.cards.last_mut().expect("a continuation follows an entry").1.extend(fields);
+        } else {
+            r.cards.push((first.trim().trim_end_matches('*').to_string(), fields));
+        }
+    }
+    r
+}
+
+/// A Nastran real: an exponent may come without `E`.
+fn real(s: &str) -> f64 {
+    let split = s.char_indices().skip(1).find(|(i, c)| (*c == '+' || *c == '-') && !s[..*i].ends_with(['E', 'e']));
+    match split {
+        Some((i, _)) => format!("{}E{}", &s[..i], &s[i..]).parse().unwrap(),
+        None => s.parse().unwrap(),
+    }
+}
+
+fn int(s: &str) -> u32 {
+    s.parse().unwrap_or_else(|_| panic!("`{s}` is not an integer"))
+}
+
+impl Read {
+    fn cards<'a>(&'a self, name: &'a str) -> impl Iterator<Item = &'a Vec<String>> + 'a {
+        self.cards.iter().filter(move |(n, _)| n == name).map(|(_, f)| f)
+    }
+
+    fn grids(&self) -> BTreeMap<u32, [f64; 3]> {
+        self.cards("GRID").map(|f| (int(&f[0]), [real(&f[2]), real(&f[3]), real(&f[4])])).collect()
+    }
+
+    /// Element connectivity in Nastran order, by label (solids and shells).
+    fn elements(&self) -> BTreeMap<u32, (String, Vec<u32>)> {
+        let mut out = BTreeMap::new();
+        for (name, f) in &self.cards {
+            let n = match name.as_str() {
+                "CTETRA" | "CPENTA" | "CHEXA" => f.len() - 2,
+                "CTRIA3" => 3,
+                "CQUAD4" => 4,
+                "CTRIA6" => 6,
+                "CQUAD8" => 8,
+                "CBEAM" | "CROD" => 2,
+                _ => continue,
+            };
+            let nodes = f[2..2 + n].iter().filter(|s| !s.is_empty()).map(|s| int(s)).collect();
+            out.insert(int(&f[0]), (name.clone(), nodes));
+        }
+        out
+    }
+
+    /// Constraint set `sid`: node → (degree of freedom → value).
+    fn constraints(&self, sid: u32) -> BTreeMap<u32, BTreeMap<u8, f64>> {
+        let mut out: BTreeMap<u32, BTreeMap<u8, f64>> = BTreeMap::new();
+        let dofs = |c: &str| c.bytes().map(|b| b - b'0').collect::<Vec<u8>>();
+        for f in self.cards("SPC1").filter(|f| int(&f[0]) == sid) {
+            for g in f[2..].iter().filter(|s| !s.is_empty()) {
+                for d in dofs(&f[1]) {
+                    assert!(out.entry(int(g)).or_default().insert(d, 0.0).is_none(), "constrained twice");
+                }
+            }
+        }
+        for f in self.cards("SPC").filter(|f| int(&f[0]) == sid) {
+            for d in dofs(&f[2]) {
+                assert!(out.entry(int(&f[1])).or_default().insert(d, real(&f[3])).is_none(), "constrained twice");
+            }
+        }
+        out
+    }
+}
+
+// ---------------------------------------------------------------------------------------------
+// Our own models.
+
+fn steel() -> Material {
+    Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9) }
+}
+
+/// A bar of `n` unit hex8 along X; node (i, y, z) = 1 + 4i + 2y + z.
+fn bar(n: u32) -> Model {
+    let mut m = Model::default();
+    m.set_units(UnitSystem::MmTonneS);
+    for i in 0..=n {
+        for (j, (y, z)) in [(0.0, 0.0), (0.0, 1.0), (1.0, 0.0), (1.0, 1.0)].into_iter().enumerate() {
+            m.set_node(1 + 4 * i + j as u32, [f64::from(i), y, z]).unwrap();
+        }
+    }
+    let g = |i: u32, y: u32, z: u32| 1 + 4 * i + 2 * y + z;
+    for e in 0..n {
+        let nodes = vec![g(e, 0, 0), g(e + 1, 0, 0), g(e + 1, 1, 0), g(e, 1, 0), g(e, 0, 1), g(e + 1, 0, 1), g(e + 1, 1, 1), g(e, 1, 1)];
+        m.set_element(e + 1, ElementKind::Hex8, nodes).unwrap();
+    }
+    m.add_set("all", Set { elements: (1..=n).collect(), ..Set::default() }).unwrap();
+    m.add_set("root", Set { faces: vec![ElementFace { element: 1, face: 6 }], ..Set::default() }).unwrap();
+    m.add_set("end", Set { faces: vec![ElementFace { element: n, face: 4 }], ..Set::default() }).unwrap();
+    m.add_material("Steel", steel()).unwrap();
+    m.add_section("bar", Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() })
+        .unwrap();
+    m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("root".into()) }).unwrap();
+    m
+}
+
+fn static_step(name: &str, bcs: &[&str], loads: &[&str]) -> Step {
+    let mut s = Step::new(name, Procedure::Static {});
+    s.boundary_conditions = bcs.iter().map(|s| s.to_string()).collect();
+    s.loads = loads.iter().map(|s| s.to_string()).collect();
+    s
+}
+
+fn cantilever() -> Model {
+    let mut m = bar(3);
+    m.add_load("tip", Load::Force { target: Target::Set("end".into()), vector: [0.0, 0.0, -100.0], distribution: Distribution::Total }).unwrap();
+    m.add_step(static_step("Bend", &["clamp"], &["tip"]), None).unwrap();
+    m
+}
+
+/// A 2 × 2 quad4 plate in z = 0, 10 × 10, clamped along x = 0, under pressure and gravity.
+fn plate() -> Model {
+    let mut m = Model::default();
+    m.set_units(UnitSystem::MmTonneS);
+    for j in 0..3u32 {
+        for i in 0..3u32 {
+            m.set_node(1 + i + 3 * j, [5.0 * f64::from(i), 5.0 * f64::from(j), 0.0]).unwrap();
+        }
+    }
+    for (e, (i, j)) in [(0, 0), (1, 0), (0, 1), (1, 1)].into_iter().enumerate() {
+        let n = 1 + i + 3 * j;
+        m.set_element(e as u32 + 1, ElementKind::Quad4, vec![n, n + 1, n + 4, n + 3]).unwrap();
+    }
+    m.add_set("all", Set { elements: vec![1, 2, 3, 4], ..Set::default() }).unwrap();
+    m.add_set("edge", Set { nodes: vec![1, 4, 7], ..Set::default() }).unwrap();
+    m.add_set("top", Set { faces: (1..=4).map(|element| ElementFace { element, face: 1 }).collect(), ..Set::default() }).unwrap();
+    m.add_material("Aluminium", Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9) }).unwrap();
+    let shell = SectionKind::Shell { thickness: 2.0, offset: 0.0, points: 5 };
+    m.add_section("skin", Section { set: "all".into(), material: "Aluminium".into(), kind: shell, formulation: Formulation::default() }).unwrap();
+    m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("edge".into()) }).unwrap();
+    m.add_load("press", 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(static_step("Sag", &["clamp"], &["press", "weight"]), None).unwrap();
+    m
+}
+
+/// Reference corners of each solid kind, slightly distorted, with straight edges (mid-side nodes
+/// at edge midpoints, by the model's own edge table).
+fn solid(kind: ElementKind) -> Vec<[f64; 3]> {
+    let corners: Vec<[f64; 3]> = match kind {
+        ElementKind::Tet4 | ElementKind::Tet10 => vec![[0.0, 0.0, 0.0], [2.0, 0.1, 0.0], [0.2, 1.5, 0.1], [0.3, 0.2, 1.7]],
+        ElementKind::Wedge6 => vec![[0.0, 0.0, 0.0], [2.0, 0.0, 0.1], [0.1, 1.5, 0.0], [0.1, 0.1, 1.2], [2.1, 0.0, 1.3], [0.0, 1.6, 1.1]],
+        _ => vec![
+            [0.0, 0.0, 0.0],
+            [2.0, 0.1, 0.0],
+            [2.2, 1.4, 0.1],
+            [-0.1, 1.5, 0.0],
+            [0.1, 0.0, 1.0],
+            [2.0, 0.2, 1.1],
+            [2.1, 1.6, 1.2],
+            [0.0, 1.4, 0.9],
+        ],
+    };
+    let mut xyz = corners.clone();
+    xyz.resize(kind.node_count(), [0.0; 3]);
+    for e in kind.edges() {
+        if let [a, b, mid] = **e {
+            xyz[mid] = [0, 1, 2].map(|i| (corners[a][i] + corners[b][i]) / 2.0);
+        }
+    }
+    xyz
+}
+
+/// One element of `kind` with its nodes labelled in reverse (so labels and positions differ), a
+/// section, and a static step.
+fn one_element(kind: ElementKind) -> Model {
+    let mut m = Model::default();
+    let xyz = solid(kind);
+    let n = xyz.len() as u32;
+    for (i, p) in xyz.iter().enumerate() {
+        m.set_node(100 - i as u32, *p).unwrap();
+    }
+    m.set_element(7, kind, (0..n).map(|i| 100 - i).collect()).unwrap();
+    m.add_set("all", Set { elements: vec![7], ..Set::default() }).unwrap();
+    m.add_set("base", Set { faces: vec![ElementFace { element: 7, face: 1 }], ..Set::default() }).unwrap();
+    m.add_material("Steel", steel()).unwrap();
+    m.add_section("s", Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() })
+        .unwrap();
+    m.add_boundary_condition("hold", BoundaryCondition::Fixed { target: Target::Set("base".into()) }).unwrap();
+    m.add_load("push", Load::Pressure { set: "base".into(), magnitude: 1.0 }).unwrap();
+    m.add_step(static_step("Load", &["hold"], &["push"]), None).unwrap();
+    m
+}
+
+/// A hex20 block whose top face is coupled to a reference point, clamped at the base, in a modes
+/// step.
+fn modes() -> Model {
+    let mut m = one_element(ElementKind::Hex20);
+    m.add_set("top", Set { faces: vec![ElementFace { element: 7, face: 2 }], ..Set::default() }).unwrap();
+    m.add_reference_point("RP", ReferencePoint { xyz: [1.0, 0.75, 2.0], label: 0 }).unwrap();
+    m.add_coupling("lid", Coupling { reference_point: "RP".into(), set: "top".into(), kind: CouplingKind::Kinematic, dofs: Dofs::all() }).unwrap();
+    m.remove_step("Load").unwrap();
+    let mut s = Step::new("Modes", Procedure::Frequency { modes: 6, min_hz: None, max_hz: Some(5000.0) });
+    s.boundary_conditions.insert("hold".into());
+    m.add_step(s, None).unwrap();
+    m
+}
+
+fn golden(name: &str, deck: &str) {
+    let path = std::path::Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/golden/nastran").join(format!("{name}.bdf"));
+    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 want = std::fs::read_to_string(&path).unwrap_or_else(|e| panic!("{}: {e} (bless with MODELCRAFT_BLESS=1)", path.display()));
+    assert_eq!(deck, want, "{name}: the deck changed (bless with MODELCRAFT_BLESS=1 after checking it)");
+}
+
+fn unsupported(m: &Model) -> String {
+    match nastran::write(m) {
+        Err(Error::Unsupported { format: "Nastran", what }) => what,
+        other => panic!("expected unsupported, got {other:?}"),
+    }
+}
+
+// ---------------------------------------------------------------------------------------------
+
+#[test]
+fn golden_decks() {
+    for (name, model) in [("cantilever", cantilever()), ("plate", plate()), ("modes", modes())] {
+        let deck = nastran::write(&model).unwrap();
+        assert_eq!(nastran::write(&model).unwrap(), deck, "deterministic");
+        golden(name, &deck);
+    }
+}
+
+#[test]
+fn cantilever_reads_back() {
+    let m = cantilever();
+    let r = read(&nastran::write(&m).unwrap());
+    assert_eq!(r.sol, 101);
+    assert_eq!(r.subcases.len(), 1);
+    let case = &r.subcases[0];
+    assert_eq!(case["SUBTITLE"], "Bend");
+    for request in ["DISPLACEMENT", "SPCFORCES", "STRESS", "STRAIN"] {
+        assert_eq!(case[request], "ALL");
+    }
+    // Grids are exact.
+    let grids = r.grids();
+    assert_eq!(grids.len(), m.nodes().len());
+    for (label, node) in m.nodes() {
+        assert_eq!(grids[label], node.xyz);
+    }
+    // Elements in Nastran order map back to the model's.
+    let elements = r.elements();
+    for (label, el) in m.elements() {
+        let (card, nodes) = &elements[label];
+        assert_eq!(card, "CHEXA");
+        assert_eq!(*nodes, el.nodes);
+    }
+    // The clamp holds every degree of freedom of the root face's nodes.
+    let spc = r.constraints(int(&case["SPC"]));
+    assert_eq!(spc.keys().copied().collect::<Vec<_>>(), vec![1, 2, 3, 4]);
+    assert!(spc.values().all(|d| d.len() == 6 && d.values().all(|v| *v == 0.0)));
+    // The tip force totals -100 in Z, a quarter on each end node.
+    let sid = int(&case["LOAD"]);
+    let mut total = [0.0; 3];
+    for f in r.cards("FORCE").filter(|f| int(&f[0]) == sid) {
+        assert!((13..=16).contains(&int(&f[1])));
+        for k in 0..3 {
+            total[k] += real(&f[3]) * real(&f[4 + k]);
+        }
+    }
+    assert_eq!(total, [0.0, 0.0, -100.0]);
+    // Two-point full integration by default.
+    let psolid = r.cards("PSOLID").next().unwrap();
+    assert_eq!((psolid[3].as_str(), psolid[5].as_str()), ("TWO", "FULL"));
+}
+
+#[test]
+fn plate_reads_back() {
+    let r = read(&nastran::write(&plate()).unwrap());
+    let case = &r.subcases[0];
+    // Pressure and gravity combine through a LOAD entry.
+    let load = r.cards("LOAD").next().unwrap();
+    assert_eq!(int(&load[0]), int(&case["LOAD"]));
+    let parts: Vec<u32> = load[2..].chunks(2).filter(|c| !c[1].is_empty()).map(|c| int(&c[1])).collect();
+    assert_eq!(parts.len(), 2);
+    let grav = r.cards("GRAV").next().unwrap();
+    assert!(parts.contains(&int(&grav[0])));
+    assert_eq!([real(&grav[2]) * real(&grav[3]), real(&grav[4]), real(&grav[5])], [0.0, 0.0, -9810.0]);
+    // Pressure on the side the normal points to pushes against the normal.
+    let pload: Vec<_> = r.cards("PLOAD4").collect();
+    assert_eq!(pload.len(), 4);
+    assert!(pload.iter().all(|f| real(&f[2]) == -0.01 && parts.contains(&int(&f[0]))));
+    let pshell = r.cards("PSHELL").next().unwrap();
+    assert_eq!(real(&pshell[2]), 2.0);
+    assert_eq!(r.elements().values().filter(|(c, _)| c == "CQUAD4").count(), 4);
+    assert_eq!(r.constraints(int(&case["SPC"])).keys().copied().collect::<Vec<_>>(), vec![1, 4, 7]);
+}
+
+#[test]
+fn modes_reads_back() {
+    let m = modes();
+    let r = read(&nastran::write(&m).unwrap());
+    assert_eq!(r.sol, 103);
+    let case = &r.subcases[0];
+    assert!(!case.contains_key("LOAD"));
+    assert_eq!(case["DISPLACEMENT"], "ALL");
+    let eigrl = r.cards("EIGRL").next().unwrap();
+    assert_eq!(int(&eigrl[0]), int(&case["METHOD"]));
+    assert_eq!((eigrl[1].as_str(), real(&eigrl[2]), int(&eigrl[3])), ("", 5000.0, 6));
+    // The reference point is a grid; the coupling a rigid element on the top face's eight nodes.
+    let rp = m.reference_point("RP").unwrap().label;
+    assert!(r.grids().contains_key(&rp));
+    let rbe2 = r.cards("RBE2").next().unwrap();
+    assert_eq!((int(&rbe2[0]), int(&rbe2[1]), int(&rbe2[2])), (8, rp, 123456));
+    let mut dependent: Vec<u32> = rbe2[3..].iter().filter(|s| !s.is_empty()).map(|s| int(s)).collect();
+    dependent.sort_unstable();
+    assert_eq!(dependent, m.resolve_set("top").unwrap().nodes);
+    // The CHEXA's mid-side blocks are swapped and map back.
+    let (_, nodes) = &r.elements()[&7];
+    let el = m.element(7).unwrap();
+    let back: Vec<u32> = (0..20).map(|i| nodes[nastran::order(ElementKind::Hex20).iter().position(|&p| p == i).unwrap()]).collect();
+    assert_eq!(back, el.nodes);
+    assert_ne!(*nodes, el.nodes);
+    assert!(r.cards("PARAM").any(|f| f[0] == "COUPMASS"));
+}
+
+// ---------------------------------------------------------------------------------------------
+// Node order against Nastran's own conventions.
+
+/// Nastran's natural coordinates of each grid, from its corner layout and its edge table (which
+/// grids sit between which corners), written independently of the writer's order table.
+fn nastran_natural(card: &str, n: usize) -> Vec<[f64; 3]> {
+    let (corners, edges): (Vec<[f64; 3]>, Vec<(usize, usize)>) = match card {
+        "CHEXA" => (
+            vec![
+                [-1.0, -1.0, -1.0],
+                [1.0, -1.0, -1.0],
+                [1.0, 1.0, -1.0],
+                [-1.0, 1.0, -1.0],
+                [-1.0, -1.0, 1.0],
+                [1.0, -1.0, 1.0],
+                [1.0, 1.0, 1.0],
+                [-1.0, 1.0, 1.0],
+            ],
+            // G9–G12 bottom, G13–G16 vertical, G17–G20 top.
+            vec![(1, 2), (2, 3), (3, 4), (4, 1), (1, 5), (2, 6), (3, 7), (4, 8), (5, 6), (6, 7), (7, 8), (8, 5)],
+        ),
+        // (ξ, η) in the triangle, ζ through the thickness; G7–G9 one triangle, G10–G12 vertical.
+        "CPENTA" => (
+            vec![[0.0, 0.0, -1.0], [1.0, 0.0, -1.0], [0.0, 1.0, -1.0], [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0]],
+            vec![(1, 2), (2, 3), (3, 1), (1, 4), (2, 5), (3, 6), (4, 5), (5, 6), (6, 4)],
+        ),
+        // G5–G10 on edges 1-2, 2-3, 3-1, 1-4, 2-4, 3-4.
+        _ => (vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], vec![(1, 2), (2, 3), (3, 1), (1, 4), (2, 4), (3, 4)]),
+    };
+    let mut out = corners.clone();
+    for (a, b) in edges.into_iter().take(n - corners.len()) {
+        out.push([0, 1, 2].map(|i| (corners[a - 1][i] + corners[b - 1][i]) / 2.0));
+    }
+    out
+}
+
+/// Linear (corner) shape functions of each Nastran solid.
+fn corner_shape(card: &str, p: [f64; 3]) -> Vec<f64> {
+    let [x, y, z] = p;
+    match card {
+        "CHEXA" => nastran_natural("CHEXA", 8).iter().map(|c| (1.0 + x * c[0]) * (1.0 + y * c[1]) * (1.0 + z * c[2]) / 8.0).collect(),
+        "CPENTA" => {
+            let l = [1.0 - x - y, x, y];
+            [(1.0 - z) / 2.0, (1.0 + z) / 2.0].iter().flat_map(|h| l.map(|li| li * h)).collect()
+        }
+        _ => vec![1.0 - x - y - z, x, y, z],
+    }
+}
+
+/// The Jacobian determinant of the corner map at natural point `p`, by central differences.
+fn jacobian(card: &str, corners: &[[f64; 3]], p: [f64; 3]) -> f64 {
+    let h = 1e-6;
+    let at = |q: [f64; 3]| -> [f64; 3] {
+        let n = corner_shape(card, q);
+        [0, 1, 2].map(|k| n.iter().zip(corners).map(|(w, c)| w * c[k]).sum())
+    };
+    let d: Vec<[f64; 3]> = (0..3)
+        .map(|axis| {
+            let (mut a, mut b) = (p, p);
+            a[axis] += h;
+            b[axis] -= h;
+            let (fa, fb) = (at(a), at(b));
+            [0, 1, 2].map(|k| (fa[k] - fb[k]) / (2.0 * h))
+        })
+        .collect();
+    d[0][0] * (d[1][1] * d[2][2] - d[1][2] * d[2][1]) - d[0][1] * (d[1][0] * d[2][2] - d[1][2] * d[2][0])
+        + d[0][2] * (d[1][0] * d[2][1] - d[1][1] * d[2][0])
+}
+
+/// Check a written solid against Nastran's conventions: positive Jacobian at points near every
+/// corner and at the centre, and every mid-side grid at the midpoint of its Nastran edge.
+fn check_nastran_solid(card: &str, xyz: &[[f64; 3]]) -> Result<(), String> {
+    let natural = nastran_natural(card, xyz.len());
+    let ncorner = corner_shape(card, [0.0; 3]).len();
+    let corners = &xyz[..ncorner];
+    let centre: [f64; 3] = [0, 1, 2].map(|k| natural[..ncorner].iter().map(|c| c[k]).sum::<f64>() / ncorner as f64);
+    for c in natural[..ncorner].iter().chain(std::iter::once(&centre)) {
+        let p = [0, 1, 2].map(|k| centre[k] + 0.9 * (c[k] - centre[k]));
+        let j = jacobian(card, corners, p);
+        if !(j > 0.0) {
+            return Err(format!("{card}: Jacobian {j} at {p:?}"));
+        }
+    }
+    for (i, nat) in natural.iter().enumerate().skip(ncorner) {
+        let n = corner_shape(card, *nat);
+        let want: [f64; 3] = [0, 1, 2].map(|k| n.iter().zip(corners).map(|(w, c)| w * c[k]).sum());
+        if (0..3).any(|k| (want[k] - xyz[i][k]).abs() > 1e-12) {
+            return Err(format!("{card}: grid G{} is not on its edge", i + 1));
+        }
+    }
+    Ok(())
+}
+
+#[test]
+fn solids_keep_nastran_node_order_and_positive_jacobians() {
+    for kind in [ElementKind::Tet4, ElementKind::Tet10, ElementKind::Wedge6, ElementKind::Hex8, ElementKind::Hex20] {
+        let m = one_element(kind);
+        let r = read(&nastran::write(&m).unwrap());
+        let grids = r.grids();
+        let (card, nodes) = &r.elements()[&7];
+        assert_eq!(card, nastran::element_card(kind));
+        let xyz: Vec<[f64; 3]> = nodes.iter().map(|n| grids[n]).collect();
+        check_nastran_solid(card, &xyz).unwrap_or_else(|e| panic!("{kind:?}: {e}"));
+        // Mapped back, the connectivity is the model's, and so is its positive Jacobian.
+        let order = nastran::order(kind);
+        let back: Vec<u32> = (0..nodes.len()).map(|i| nodes[order.iter().position(|&p| p == i).unwrap()]).collect();
+        assert_eq!(back, m.element(7).unwrap().nodes, "{kind:?}");
+
+        // The face pressure names a corner on the face, and the corner diagonally opposite (hex),
+        // the corner off the face (tet), or nothing (a wedge's triangle).
+        let p = r.cards("PLOAD4").next().unwrap();
+        assert_eq!(real(&p[2]), 1.0);
+        let face: Vec<u32> = m.face_nodes(ElementFace { element: 7, face: 1 }).unwrap();
+        assert!(face.contains(&int(&p[6])));
+        match kind {
+            ElementKind::Tet4 | ElementKind::Tet10 => assert!(!face.contains(&int(&p[7])) && nodes[..4].contains(&int(&p[7]))),
+            ElementKind::Wedge6 => assert_eq!(p.get(7).map(String::as_str).unwrap_or(""), ""),
+            _ => assert_eq!(int(&p[7]), face[2]),
+        }
+    }
+    // Without the reordering, a hex20 breaks Nastran's convention: the check catches it.
+    let m = one_element(ElementKind::Hex20);
+    let model_order: Vec<[f64; 3]> = m.element(7).unwrap().nodes.iter().map(|n| m.node(*n).unwrap().xyz).collect();
+    assert!(check_nastran_solid("CHEXA", &model_order).unwrap_err().contains("G13"));
+}
+
+// ---------------------------------------------------------------------------------------------
+// The rest of the model.
+
+#[test]
+fn beams_rods_orthotropy_couplings_and_prescribed_motion() {
+    let mut m = cantilever();
+    // An orthotropic solid material.
+    let ortho = Material {
+        elastic: Some(Elastic::EngineeringConstants { young: [140e3, 10e3, 12e3], poisson: [0.3, 0.28, 0.4], shear: [5e3, 4.5e3, 3.8e3] }),
+        density: Some(1.6e-9),
+    };
+    m.add_material("Carbon", ortho.clone()).unwrap();
+    m.replace_section(
+        "bar",
+        Section { set: "all".into(), material: "Carbon".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() },
+    )
+    .unwrap();
+    // A beam and a rod hanging off the tip, and a shell lying in z = 0 with the same material.
+    m.set_node(50, [5.0, 0.0, 0.0]).unwrap();
+    m.set_node(51, [4.0, 1.0, 0.0]).unwrap();
+    m.set_element(10, ElementKind::Beam2, vec![13, 50]).unwrap();
+    m.set_element(11, ElementKind::Rod2, vec![15, 51]).unwrap();
+    m.set_element(12, ElementKind::Tri3, vec![13, 50, 51]).unwrap();
+    m.add_set("beam", Set { elements: vec![10], ..Set::default() }).unwrap();
+    m.add_set("rod", Set { elements: vec![11], ..Set::default() }).unwrap();
+    m.add_set("web", Set { elements: vec![12], ..Set::default() }).unwrap();
+    let i = BeamProfile::I { height: 10.0, width: 6.0, web: 1.0, flange: 1.5 };
+    m.add_section(
+        "beam",
+        Section {
+            set: "beam".into(),
+            material: "Steel".into(),
+            kind: SectionKind::Beam { profile: i, orientation: [0.0, 0.0, 1.0] },
+            formulation: Formulation::default(),
+        },
+    )
+    .unwrap();
+    m.add_section(
+        "rod",
+        Section { set: "rod".into(), material: "Steel".into(), kind: SectionKind::Rod { area: 3.5 }, formulation: Formulation::default() },
+    )
+    .unwrap();
+    m.add_section(
+        "web",
+        Section {
+            set: "web".into(),
+            material: "Carbon".into(),
+            kind: SectionKind::Shell { thickness: 1.0, offset: 0.5, points: 3 },
+            formulation: Formulation::default(),
+        },
+    )
+    .unwrap();
+    // A distributing coupling from a reference point to the end face, loaded with a moment.
+    m.add_reference_point("Hub", ReferencePoint { xyz: [3.0, 0.5, 0.5], label: 0 }).unwrap();
+    m.add_coupling("spider", Coupling { reference_point: "Hub".into(), set: "end".into(), kind: CouplingKind::Distributing, dofs: Dofs::all() })
+        .unwrap();
+    m.add_load("twist", Load::Moment { target: Target::ReferencePoint("Hub".into()), vector: [1000.0, 0.0, 0.0], distribution: Distribution::Total })
+        .unwrap();
+    // Prescribed motion of the root in a second step.
+    m.add_boundary_condition(
+        "lift",
+        BoundaryCondition::Displacement { target: Target::Set("root".into()), values: [Some(0.0), Some(0.0), Some(0.25), None, None, None] },
+    )
+    .unwrap();
+    m.add_step(static_step("Lift", &["lift"], &["twist"]), None).unwrap();
+
+    let deck = nastran::write(&m).unwrap();
+    let r = read(&deck);
+    assert_eq!(r.subcases.len(), 2);
+
+    // MAT9 carries the stiffness of the engineering constants; MAT8 the shell's constants.
+    let mat9 = r.cards("MAT9").next().unwrap();
+    let c = nastran::orthotropic_stiffness([140e3, 10e3, 12e3], [0.3, 0.28, 0.4], [5e3, 4.5e3, 3.8e3]).unwrap();
+    let close = |s: &str, v: f64| (real(s) - v).abs() <= 1e-13 * v.abs();
+    assert!(close(&mat9[1], c[0][0]) && close(&mat9[7], c[1][0]) && close(&mat9[12], c[1][2]));
+    assert_eq!((real(&mat9[16]), real(&mat9[19]), real(&mat9[21])), (5e3, 3.8e3, 4.5e3));
+    assert_eq!(real(&mat9[22]), 1.6e-9);
+    let mat8 = r.cards("MAT8").next().unwrap();
+    assert_eq!(mat8[1..7].iter().map(|s| real(s)).collect::<Vec<_>>(), vec![140e3, 10e3, 0.3, 5e3, 4.5e3, 3.8e3]);
+    // The shell's material axes are the basic ones, and its offset puts the mid-surface below the nodes.
+    let tri = r.cards("CTRIA3").next().unwrap();
+    assert_eq!((int(&tri[5]), real(&tri[6])), (1, -0.5));
+    assert!(r.cards("CORD2R").next().is_some());
+
+    // The beam: its I section's flanges lie along the orientation (Nastran's element z axis).
+    let pbeaml = r.cards("PBEAML").next().unwrap();
+    assert_eq!(pbeaml[3], "I");
+    assert_eq!(pbeaml[8..14].iter().map(|s| real(s)).collect::<Vec<_>>(), vec![10.0, 6.0, 6.0, 1.0, 1.5, 1.5]);
+    let cbeam = r.cards("CBEAM").next().unwrap();
+    let v = [real(&cbeam[4]), real(&cbeam[5]), real(&cbeam[6])];
+    let axis = [2.0, 0.0, 0.0]; // node 13 (3, 0, 0) to node 50 (5, 0, 0)
+    let z = [axis[1] * v[2] - axis[2] * v[1], axis[2] * v[0] - axis[0] * v[2], axis[0] * v[1] - axis[1] * v[0]];
+    assert!(z[0] == 0.0 && z[1] == 0.0 && z[2] != 0.0, "element z along the orientation: {z:?}");
+    let prod = r.cards("PROD").next().unwrap();
+    assert_eq!(real(&prod[2]), 3.5);
+
+    // The distributing coupling: the hub depends on the end face's translations, equally weighted.
+    let rbe3 = r.cards("RBE3").next().unwrap();
+    let hub = m.reference_point("Hub").unwrap().label;
+    assert_eq!((rbe3[1].as_str(), int(&rbe3[2]), int(&rbe3[3]), real(&rbe3[4]), int(&rbe3[5])), ("", hub, 123456, 1.0, 123));
+    let moment = r.cards("MOMENT").next().unwrap();
+    assert_eq!((int(&moment[1]), real(&moment[4])), (hub, 1000.0));
+
+    // Step 2 prescribes Z = 0.25 at the root and holds X and Y.
+    let spc = r.constraints(int(&r.subcases[1]["SPC"]));
+    assert_eq!(spc.len(), 4);
+    for dofs in spc.values() {
+        assert_eq!(dofs.iter().map(|(d, v)| (*d, *v)).collect::<Vec<_>>(), vec![(1, 0.0), (2, 0.0), (3, 0.25)]);
+    }
+    assert_ne!(r.subcases[0]["SPC"], r.subcases[1]["SPC"]);
+}
+
+#[test]
+fn steps_are_chosen_and_mixed_solutions_are_refused() {
+    let mut m = cantilever();
+    let mut s = Step::new("Modes", Procedure::Frequency { modes: 4, min_hz: Some(1.0), max_hz: None });
+    s.boundary_conditions.insert("clamp".into());
+    m.add_step(s, None).unwrap();
+    assert!(unsupported(&m).contains("separate decks"));
+
+    let only = |names: &[&str]| names.iter().map(|s| s.to_string()).collect::<Vec<_>>();
+    let r = read(&nastran::write(&with_steps(&m, &only(&["Modes"])).unwrap()).unwrap());
+    assert_eq!((r.sol, r.subcases.len()), (103, 1));
+    let eigrl = r.cards("EIGRL").next().unwrap();
+    assert_eq!((real(&eigrl[1]), eigrl[2].as_str(), int(&eigrl[3])), (1.0, "", 4));
+    // A modes step holds the clamp but carries no loads.
+    assert!(r.cards("FORCE").next().is_none());
+    assert_eq!(read(&nastran::write(&with_steps(&m, &only(&["Bend"])).unwrap()).unwrap()).sol, 101);
+    let nastran = modelcraft_decks::format("nastran").unwrap();
+    assert_eq!(nastran.write_steps(&m, Some(&only(&["Bend"]))).unwrap().text, nastran::write(&with_steps(&m, &only(&["Bend"])).unwrap()).unwrap());
+
+    for bad in [only(&[]), only(&["Nope"]), only(&[&"x".repeat(100_000)]), only(&["", "\u{0}"])] {
+        assert!(matches!(with_steps(&m, &bad), Err(Error::Model(_))), "{bad:?}");
+    }
+}
+
+#[test]
+fn what_nastran_cannot_express_is_named() {
+    let solid = |f: Formulation| Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation: f };
+    let mut m = cantilever();
+    m.replace_section("bar", solid(Formulation { hybrid: true, ..Formulation::default() })).unwrap();
+    assert!(unsupported(&m).contains("hybrid"));
+    m.replace_section("bar", solid(Formulation { integration: Integration::Reduced, ..Formulation::default() })).unwrap();
+    assert!(unsupported(&m).contains("reduced integration of hex8"));
+    m.replace_section("bar", solid(Formulation { incompatible_modes: true, ..Formulation::default() })).unwrap();
+    let r = read(&nastran::write(&m).unwrap());
+    assert_eq!(r.cards("PSOLID").next().unwrap()[3], "BUBBLE");
+
+    // Several problems are reported together.
+    let mut m = cantilever();
+    m.add_load("heavy", Load::BodyForce { set: "all".into(), vector: [0.0, 0.0, -1.0] }).unwrap();
+    m.add_set("first", Set { elements: vec![1], ..Set::default() }).unwrap();
+    m.add_load("partial", Load::Gravity { acceleration: [0.0, 0.0, -9810.0], set: Some("first".into()) }).unwrap();
+    m.add_reference_point("RP", ReferencePoint { xyz: [3.0, 0.5, 0.5], label: 0 }).unwrap();
+    m.add_coupling("grip", Coupling { reference_point: "RP".into(), set: "root".into(), kind: CouplingKind::Kinematic, dofs: Dofs::all() }).unwrap();
+    m.replace_step(static_step("Bend", &["clamp"], &["tip", "heavy", "partial"])).unwrap();
+    let what = unsupported(&m);
+    for needle in ["body force", "gravity on part", "coupling `grip`"] {
+        assert!(what.contains(needle), "{needle}: {what}");
+    }
+
+    // Orthotropic beams, tilted orthotropic shells, one-point shells.
+    let mut m = plate();
+    let ortho = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] };
+    m.replace_material("Aluminium", Material { elastic: Some(ortho), density: Some(1.6e-9) }).unwrap();
+    m.set_node(9, [10.0, 10.0, 0.5]).unwrap();
+    assert!(unsupported(&m).contains("planes of constant Z"));
+    m.set_node(9, [10.0, 10.0, 0.0]).unwrap();
+    nastran::write(&m).unwrap();
+    m.replace_section(
+        "skin",
+        Section {
+            set: "all".into(),
+            material: "Aluminium".into(),
+            kind: SectionKind::Shell { thickness: 2.0, offset: 0.0, points: 1 },
+            formulation: Formulation::default(),
+        },
+    )
+    .unwrap();
+    assert!(unsupported(&m).contains("one integration point"));
+    let mut m = plate();
+    m.set_element(9, ElementKind::Rod2, vec![1, 9]).unwrap();
+    m.add_set("bars", Set { elements: vec![9], ..Set::default() }).unwrap();
+    m.add_material("Fibre", Material { elastic: Some(ortho), density: Some(1.6e-9) }).unwrap();
+    m.add_section(
+        "bars",
+        Section { set: "bars".into(), material: "Fibre".into(), kind: SectionKind::Rod { area: 1.0 }, formulation: Formulation::default() },
+    )
+    .unwrap();
+    assert!(unsupported(&m).contains("beams or rods"));
+
+    // A model that does not resolve is the model's error.
+    assert!(matches!(nastran::write(&Model::default()), Err(Error::Model(_))));
+}
+
+#[test]
+fn conflicting_values_and_ids_at_the_top_of_the_range() {
+    let mut m = cantilever();
+    m.add_boundary_condition(
+        "raise",
+        BoundaryCondition::Displacement { target: Target::Set("root".into()), values: [None, None, Some(1.0), None, None, None] },
+    )
+    .unwrap();
+    m.replace_step(static_step("Bend", &["clamp", "raise"], &["tip"])).unwrap();
+    assert!(unsupported(&m).contains("different values"));
+
+    // An element at the top of the label range leaves no id for a coupling's rigid element.
+    let mut m = cantilever();
+    m.set_node(99_999_998, [9.0, 9.0, 9.0]).unwrap();
+    m.set_element(99_999_999, ElementKind::Rod2, vec![1, 99_999_998]).unwrap();
+    m.replace_set("all", Set { elements: vec![1, 2, 3], ..Set::default() }).unwrap();
+    m.add_set("bar2", Set { elements: vec![99_999_999], ..Set::default() }).unwrap();
+    m.add_section(
+        "bar2",
+        Section { set: "bar2".into(), material: "Steel".into(), kind: SectionKind::Rod { area: 1.0 }, formulation: Formulation::default() },
+    )
+    .unwrap();
+    let deck = nastran::write(&m).unwrap();
+    assert!(deck.contains("CROD    99999999"));
+    m.add_reference_point("RP", ReferencePoint { xyz: [3.0, 0.5, 0.5], label: 1000 }).unwrap();
+    m.add_coupling("pad", Coupling { reference_point: "RP".into(), set: "end".into(), kind: CouplingKind::Kinematic, dofs: Dofs::all() }).unwrap();
+    assert!(matches!(nastran::write(&m), Err(Error::Limit(_))));
+}

crates/engine/src/cmd/deck.rs

@@ -10,7 +10,7 @@ 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"#,
+    params: r#"{"format": "abaqus|nastran", "path": "part.inp", "steps": ["Load"]} — writes the model as a solver deck (formats: modelcraft_decks::FORMATS); `steps` (optional) writes only those steps, as a Nastran run needs (static or frequency steps, not both); fails naming every error of the model check and everything the format cannot express"#,
     kind: CommandKind::Session,
     run: write,
 };
@@ -23,7 +23,8 @@ fn write(s: &mut Session, p: &Value) -> Result<Value> {
         EngineError::Params(format!("unknown deck format `{shown}`; known: {}", known.join(", ")))
     })?;
     let path = super::path_param(p, "path")?;
-    let deck = format.write(&s.model)?;
+    let steps: Option<Vec<String>> = super::opt(p, "steps")?;
+    let deck = format.write_steps(&s.model, steps.as_deref())?;
     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}))
 }
@@ -60,6 +61,20 @@ mod tests {
         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(_))));
+
+        // Nastran: one run is static or frequency steps, so the bar's two steps need two decks.
+        let bdf = dir.join("bar.bdf").to_string_lossy().into_owned();
+        let e = s.execute("deck.write", &json!({"format": "nastran", "path": bdf})).unwrap_err().to_string();
+        assert!(e.contains("not supported by Nastran") && e.contains("separate decks"), "{e}");
+        for (step, sol) in [("Bend", "SOL 101"), ("Modes", "SOL 103")] {
+            s.execute("deck.write", &json!({"format": "nastran", "path": bdf, "steps": [step]})).unwrap();
+            let deck = std::fs::read_to_string(&bdf).unwrap();
+            assert!(deck.contains(sol) && deck.contains("RBE2") && deck.ends_with("ENDDATA\n"), "{deck}");
+        }
+        assert_eq!(s.journal.len(), journal);
+        for steps in [json!([]), json!(["Nope"]), json!("Bend"), json!([1, 2]), json!([["Bend"]])] {
+            assert!(s.execute("deck.write", &json!({"format": "nastran", "path": bdf, "steps": steps})).is_err(), "{steps}");
+        }
         let _ = std::fs::remove_dir_all(&dir);
     }
 }

crates/engine/src/lib.rs

@@ -239,6 +239,11 @@ mod tests {
             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!({"format": "nastran", "path": "/no/such/dir/part.bdf"}),
+            json!({"format": "nastran", "path": "/no/such/dir/part.bdf", "steps": ["Bend"]}),
+            json!({"format": "nastran", "path": "/no/such/dir/part.bdf", "steps": [null, 3, "Modes"]}),
+            json!({"format": "nastran", "path": "/no/such/dir/part.bdf", "steps": vec!["x".repeat(1000); 2000]}),
+            json!({"format": "nastran", "path": "/no/such/dir/part.bdf", "steps": {}}),
             json!({"path": 7}),
             json!({"body": "Body1", "tolerance": 0.05}),
             json!({"body": null, "tolerance": "fine"}),

crates/model/src/bcs.rs

@@ -32,6 +32,16 @@ impl Dofs {
         self.0 == 0
     }
 
+    /// The degrees of freedom as bits: bit `d - 1` for degree of freedom `d`.
+    pub fn bits(self) -> u8 {
+        self.0
+    }
+
+    /// The degrees of freedom of bits (bits above the sixth are ignored).
+    pub fn from_bits(bits: u8) -> Dofs {
+        Dofs(bits & 0b111_111)
+    }
+
     /// The degrees of freedom in order.
     pub fn list(self) -> Vec<u8> {
         (1..=6).filter(|d| self.contains(*d)).collect()

crates/model/src/mesh.rs

@@ -76,6 +76,22 @@ impl ElementKind {
         }
     }
 
+    /// Corner nodes: the first ones in the node list (the rest are mid-side nodes).
+    pub fn corner_count(self) -> usize {
+        match self {
+            ElementKind::Rod2 | ElementKind::Beam2 => 2,
+            ElementKind::Tri3 | ElementKind::Tri6 => 3,
+            ElementKind::Quad4 | ElementKind::Quad8 | ElementKind::Tet4 | ElementKind::Tet10 => 4,
+            ElementKind::Wedge6 => 6,
+            ElementKind::Hex8 | ElementKind::Hex20 => 8,
+        }
+    }
+
+    /// Whether the element has mid-side nodes.
+    pub fn is_quadratic(self) -> bool {
+        self.node_count() > self.corner_count()
+    }
+
     /// Topological dimension: 1 for line elements, 2 for shells, 3 for solids.
     pub fn dimension(self) -> u8 {
         match self {
@@ -387,6 +403,10 @@ mod tests {
         assert_eq!(ElementKind::Tet4.faces().len(), 4);
         assert_eq!(ElementKind::Wedge6.faces().len(), 5);
         assert_eq!(ElementKind::Hex20.edges().len(), 12);
+        for kind in ElementKind::ALL {
+            let mid_sides = kind.edges().iter().filter(|e| e.len() == 3).count();
+            assert_eq!(kind.corner_count() + mid_sides, kind.node_count(), "{kind:?}");
+        }
     }
 
     #[test]

docs/parity.md

@@ -57,7 +57,7 @@ products are run, no documentation text is copied). Status: **done**, **started*
 | Interaction | Contact (general and pairs), friction, automatic pair detection | M7 | | planned |
 | Interaction | Bolt pretension, springs, point masses, connectors | M7 | | 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 | Nastran bulk data writer (`deck.write`, SOL 101 or 103 per deck; reports what Nastran cannot express) | M3 | | done |
 | Decks | LS-DYNA keyword writer (implicit; explicit in M8) | M3 | | planned |
 | Decks | Readers keeping unknown keywords | M9 | | planned |
 | Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | planned |
@@ -81,3 +81,10 @@ products are run, no documentation text is copied). Status: **done**, **started*
   `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.
+- Nastran (`modelcraft_decks::nastran`): entries with a real use large (16-column) fields, the rest
+  small fields; a real is exact when its shortest round-trip form fits 16 columns, else rounded to
+  at least ten significant digits. Linear solids map full integration to `PSOLID` `TWO`/`THREE` +
+  `FULL`, hex20 reduced to `TWO` + `FULL`, hex8 incompatible modes to `BUBBLE`; one-point reduced
+  linear solids, hybrid solids, reduced-integration or one-point shells, body forces, gravity on part
+  of the model, orthotropic beams and rods, and orthotropic shells off planes of constant Z are
+  reported as unsupported. Static and frequency steps need separate decks (`deck.write` `steps`).