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

//! Abaqus writer tests: golden decks of our own small models, a re-read of the written mesh and
//! sets by a minimal test-only reader, element orientation (Jacobian sign) from the written
//! connectivity, number round trips, line limits and everything reported as unsupported.

use std::collections::{BTreeMap, BTreeSet};
use std::path::PathBuf;

use modelcraft_decks::{Error, abaqus};
use modelcraft_model::{
    BeamProfile, BoundaryCondition, Coupling, CouplingKind, Distribution, Dofs, Elastic, ElementFace, ElementKind, Formulation, Integration, Load,
    Material, Model, OutputVariable, Procedure, ReferencePoint, Section, SectionKind, Set, Step, Symmetry, Target, UnitSystem,
};

// ---------------------------------------------------------------------------------------------
// Models

fn steel() -> Material {
    Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), fluid: None }
}

fn faces(list: &[(u32, u8)]) -> Set {
    Set { faces: list.iter().map(|&(element, face)| ElementFace { element, face }).collect(), ..Set::default() }
}

fn elements(list: impl IntoIterator<Item = u32>) -> Set {
    Set { elements: list.into_iter().collect(), ..Set::default() }
}

fn step(name: &str, procedure: Procedure, bcs: &[&str], loads: &[&str]) -> Step {
    let mut s = Step::new(name, procedure);
    s.boundary_conditions = bcs.iter().map(|s| s.to_string()).collect();
    s.loads = loads.iter().map(|s| s.to_string()).collect();
    s
}

/// A 100 × 10 × 10 cantilever of four hex8 elements along X: node (i, j, k) at
/// (25 i, 10 j, 10 k) is labelled 1 + 4 i + 2 j + k.
fn cantilever_mesh() -> Model {
    let mut m = Model::default();
    m.set_units(UnitSystem::MmTonneS);
    let n = |i: u32, j: u32, k: u32| 1 + 4 * i + 2 * j + k;
    for i in 0..=4 {
        for j in 0..=1 {
            for k in 0..=1 {
                m.set_node(n(i, j, k), [25.0 * f64::from(i), 10.0 * f64::from(j), 10.0 * f64::from(k)]).unwrap();
            }
        }
    }
    for i in 0..4 {
        let nodes = vec![n(i, 0, 0), n(i + 1, 0, 0), n(i + 1, 1, 0), n(i, 1, 0), n(i, 0, 1), n(i + 1, 0, 1), n(i + 1, 1, 1), n(i, 1, 1)];
        m.set_element(i + 1, ElementKind::Hex8, nodes).unwrap();
    }
    m.add_set("all", elements(1..=4)).unwrap();
    // Face 6 of element 1 lies on x = 0, face 4 of element 4 on x = 100.
    m.add_set("root", faces(&[(1, 6)])).unwrap();
    m.add_set("tip", faces(&[(4, 4)])).unwrap();
    m.add_material("Steel", steel()).unwrap();
    let formulation = Formulation { incompatible_modes: true, ..Formulation::default() };
    m.add_section("beam", Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation }).unwrap();
    m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("root".into()) }).unwrap();
    m
}

/// The cantilever with a total end force of 100 down.
fn cantilever() -> Model {
    let mut m = cantilever_mesh();
    m.add_load("end", Load::Force { target: Target::Set("tip".into()), vector: [0.0, 0.0, -100.0], distribution: Distribution::Total }).unwrap();
    m.add_step(step("Bend", Procedure::Static {}, &["clamp"], &["end"]), None).unwrap();
    m
}

/// A 200 × 200 plate of 2 × 2 quad4 shells in z = 0, clamped along x = 0, under pressure on its
/// top face and its own weight, with a symmetry plane along x = 100.
fn plate() -> Model {
    let mut m = Model::default();
    m.set_units(UnitSystem::MmTonneS);
    let n = |i: u32, j: u32| 1 + i + 3 * j;
    for j in 0..=2 {
        for i in 0..=2 {
            m.set_node(n(i, j), [100.0 * f64::from(i), 100.0 * f64::from(j), 0.0]).unwrap();
        }
    }
    let mut e = 1;
    for j in 0..2 {
        for i in 0..2 {
            m.set_element(e, ElementKind::Quad4, vec![n(i, j), n(i + 1, j), n(i + 1, j + 1), n(i, j + 1)]).unwrap();
            e += 1;
        }
    }
    m.add_set("skin", elements(1..=4)).unwrap();
    m.add_set("top", faces(&[(1, 1), (2, 1), (3, 1), (4, 1)])).unwrap();
    m.add_set("edge", Set { nodes: vec![1, 4, 7], ..Set::default() }).unwrap();
    m.add_set("mid", Set { nodes: vec![2, 5, 8], ..Set::default() }).unwrap();
    let aluminium = Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9), fluid: None };
    m.add_material("Aluminium", aluminium).unwrap();
    let kind = SectionKind::Shell { thickness: 2.0, offset: 0.5, points: 5 };
    let formulation = Formulation { integration: Integration::Reduced, ..Formulation::default() };
    m.add_section("plate", Section { set: "skin".into(), material: "Aluminium".into(), kind, formulation }).unwrap();
    m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("edge".into()) }).unwrap();
    m.add_boundary_condition("sym", BoundaryCondition::Symmetry { target: Target::Set("mid".into()), plane: Symmetry::X }).unwrap();
    m.add_load("wind", Load::Pressure { set: "top".into(), magnitude: 0.01 }).unwrap();
    m.add_load("weight", Load::Gravity { acceleration: [0.0, 0.0, -9810.0], set: None }).unwrap();
    m.add_step(step("Load", Procedure::Static {}, &["clamp", "sym"], &["wind", "weight"]), None).unwrap();
    m
}

/// The cantilever with a reference point coupled to its tip, three static steps that change the
/// active boundary conditions and loads, and a frequency step in between.
fn frequency() -> Model {
    let mut m = cantilever_mesh();
    m.add_reference_point("Tip", ReferencePoint { xyz: [110.0, 5.0, 5.0], label: 0 }).unwrap();
    m.add_coupling("tie", Coupling { reference_point: "Tip".into(), set: "tip".into(), kind: CouplingKind::Kinematic, dofs: Dofs::all() }).unwrap();
    let rp = || Target::ReferencePoint("Tip".into());
    m.add_load("twist", Load::Moment { target: rp(), vector: [5000.0, 0.0, 0.0], distribution: Distribution::Total }).unwrap();
    m.add_load("swell", Load::BodyForce { set: "all".into(), vector: [0.0, 1e-3, 0.0] }).unwrap();
    m.add_boundary_condition("lift", BoundaryCondition::Displacement { target: rp(), values: [None, None, Some(0.5), None, None, None] }).unwrap();
    m.add_step(step("Twist", Procedure::Static {}, &["clamp"], &["twist"]), None).unwrap();
    m.add_step(step("Modes", Procedure::Frequency { modes: 6, min_hz: Some(1.0), max_hz: None }, &["clamp", "lift"], &[]), None).unwrap();
    m.add_step(step("Lift", Procedure::Static {}, &["clamp", "lift"], &["twist", "swell"]), None).unwrap();
    let mut release = step("Release", Procedure::Static {}, &["clamp"], &[]);
    release.output = [OutputVariable::Displacement].into_iter().collect();
    m.add_step(release, None).unwrap();
    m
}

// ---------------------------------------------------------------------------------------------
// A minimal reader of what the writer emits (full deck readers are a later milestone)

#[derive(Default, Debug)]
struct Parsed {
    nodes: BTreeMap<u32, [f64; 3]>,
    elements: BTreeMap<u32, (String, Vec<u32>)>,
    nsets: BTreeMap<String, Vec<u32>>,
    elsets: BTreeMap<String, Vec<u32>>,
    surfaces: BTreeMap<String, Vec<(u32, String)>>,
}

fn param(keyword: &str, name: &str) -> Option<String> {
    keyword.split(',').skip(1).find_map(|p| {
        let (k, v) = p.split_once('=')?;
        (k.trim().eq_ignore_ascii_case(name)).then(|| v.trim().to_string())
    })
}

fn read(deck: &str) -> Parsed {
    let mut p = Parsed::default();
    let mut keyword = String::new();
    let mut pending: Vec<String> = Vec::new();
    for line in deck.lines() {
        if line.starts_with("**") {
            continue;
        }
        if line.starts_with('*') {
            keyword = line.to_string();
            continue;
        }
        let head = keyword.split(',').next().unwrap_or_default().trim().to_ascii_uppercase();
        let fields: Vec<String> = line.split(',').map(|f| f.trim().to_string()).collect();
        let label = |s: &str| s.parse::<u32>().unwrap();
        match head.as_str() {
            "*NODE" => {
                let xyz = [1, 2, 3].map(|i| fields[i].parse::<f64>().unwrap());
                p.nodes.insert(label(&fields[0]), xyz);
            }
            "*ELEMENT" => {
                pending.extend(fields.iter().filter(|f| !f.is_empty()).cloned());
                if !line.trim_end().ends_with(',') {
                    let ids: Vec<u32> = pending.drain(..).map(|f| label(&f)).collect();
                    p.elements.insert(ids[0], (param(&keyword, "TYPE").unwrap(), ids[1..].to_vec()));
                }
            }
            "*NSET" => p.nsets.entry(param(&keyword, "NSET").unwrap()).or_default().extend(fields.iter().map(|f| label(f))),
            "*ELSET" => p.elsets.entry(param(&keyword, "ELSET").unwrap()).or_default().extend(fields.iter().map(|f| label(f))),
            "*SURFACE" if param(&keyword, "TYPE").as_deref() == Some("ELEMENT") => {
                p.surfaces.entry(param(&keyword, "NAME").unwrap()).or_default().push((label(&fields[0]), fields[1].clone()));
            }
            _ => {}
        }
    }
    p
}

// ---------------------------------------------------------------------------------------------
// Golden decks

fn golden(name: &str, deck: &str) {
    let path = PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("tests/abaqus").join(format!("{name}.inp"));
    if std::env::var_os("MODELCRAFT_BLESS").is_some() {
        std::fs::create_dir_all(path.parent().unwrap()).unwrap();
        std::fs::write(&path, deck).unwrap();
    }
    let expected = std::fs::read_to_string(&path).unwrap_or_else(|e| panic!("{}: {e} (run with MODELCRAFT_BLESS=1 to write it)", path.display()));
    assert!(expected == deck, "{name}.inp differs from the written deck (MODELCRAFT_BLESS=1 rewrites it):\n{deck}");
}

fn assert_lines_fit(deck: &str) {
    for (i, line) in deck.lines().enumerate() {
        assert!(line.len() <= abaqus::MAX_LINE, "line {} has {} characters", i + 1, line.len());
        assert!(line.is_ascii(), "line {} is not ASCII", i + 1);
    }
}

#[test]
fn golden_decks() {
    for (name, model) in [("cantilever", cantilever()), ("plate", plate()), ("frequency", frequency())] {
        let deck = abaqus::write(&model).unwrap();
        assert_lines_fit(&deck);
        golden(name, &deck);
        // Deterministic: the same model, and the model after a file round trip, give the same bytes.
        assert_eq!(abaqus::write(&model).unwrap(), deck);
        let reread: Model = serde_json::from_str(&serde_json::to_string(&model).unwrap()).unwrap();
        assert_eq!(abaqus::write(&reread).unwrap(), deck, "{name}");
    }
}

#[test]
fn the_end_force_is_shared_by_area_and_keeps_its_total() {
    let deck = abaqus::write(&cantilever()).unwrap();
    let cload: Vec<&str> = deck.lines().skip_while(|l| !l.starts_with("*CLOAD")).skip(1).take_while(|l| !l.starts_with('*')).collect();
    // Four equal corners of a flat square face.
    assert_eq!(cload, ["17, 3, -25.0", "18, 3, -25.0", "19, 3, -25.0", "20, 3, -25.0"]);
    let total: f64 = cload.iter().map(|l| l.rsplit(", ").next().unwrap().parse::<f64>().unwrap()).sum();
    assert_eq!(total, -100.0);
}

#[test]
fn the_written_mesh_and_sets_read_back() {
    for model in [cantilever(), plate(), frequency()] {
        let deck = abaqus::write(&model).unwrap();
        let p = read(&deck);
        let rps: BTreeMap<u32, [f64; 3]> = model.reference_points().values().map(|r| (r.label, r.xyz)).collect();
        let mut nodes: BTreeMap<u32, [f64; 3]> = model.nodes().iter().map(|(l, n)| (*l, n.xyz)).collect();
        nodes.extend(rps);
        assert_eq!(p.nodes, nodes);
        let connectivity: BTreeMap<u32, Vec<u32>> = p.elements.iter().map(|(l, (_, n))| (*l, n.clone())).collect();
        assert_eq!(connectivity, model.elements().iter().map(|(l, e)| (*l, e.nodes.clone())).collect());
        let r = model.resolve().unwrap();
        for (name, s) in &r.sets {
            assert_eq!(p.nsets.get(name).cloned().unwrap_or_default(), s.nodes, "nset {name}");
            assert_eq!(p.elsets.get(name).cloned().unwrap_or_default(), s.elements, "elset {name}");
            let faces: Vec<(u32, String)> =
                s.faces.iter().map(|f| (f.element, abaqus::face_id(model.element(f.element).unwrap().kind, f.face).unwrap().to_string())).collect();
            assert_eq!(p.surfaces.get(name).cloned().unwrap_or_default(), faces, "surface {name}");
        }
    }
}

// ---------------------------------------------------------------------------------------------
// Node order: the element Jacobian of the written connectivity, in Abaqus's own conventions

/// Natural coordinates of each node, by Abaqus element family.
fn natural(kind: ElementKind) -> Vec<[f64; 3]> {
    let hex = [[-1., -1., -1.], [1., -1., -1.], [1., 1., -1.], [-1., 1., -1.], [-1., -1., 1.], [1., -1., 1.], [1., 1., 1.], [-1., 1., 1.]];
    match kind {
        ElementKind::Hex8 => hex.to_vec(),
        ElementKind::Hex20 => {
            let mut v = hex.to_vec();
            // Mid-sides: the bottom face's edges, the top face's, then the vertical ones.
            v.extend([[0., -1., -1.], [1., 0., -1.], [0., 1., -1.], [-1., 0., -1.]]);
            v.extend([[0., -1., 1.], [1., 0., 1.], [0., 1., 1.], [-1., 0., 1.]]);
            v.extend([[-1., -1., 0.], [1., -1., 0.], [1., 1., 0.], [-1., 1., 0.]]);
            v
        }
        ElementKind::Tet4 | ElementKind::Tet10 => {
            let mut v = vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]];
            if kind == ElementKind::Tet10 {
                v.extend([[0.5, 0., 0.], [0.5, 0.5, 0.], [0., 0.5, 0.], [0., 0., 0.5], [0.5, 0., 0.5], [0., 0.5, 0.5]]);
            }
            v
        }
        ElementKind::Wedge6 => vec![[0., 0., -1.], [1., 0., -1.], [0., 1., -1.], [0., 0., 1.], [1., 0., 1.], [0., 1., 1.]],
        ElementKind::Quad4 => vec![[-1., -1., 0.], [1., -1., 0.], [1., 1., 0.], [-1., 1., 0.]],
        ElementKind::Quad8 => {
            vec![[-1., -1., 0.], [1., -1., 0.], [1., 1., 0.], [-1., 1., 0.], [0., -1., 0.], [1., 0., 0.], [0., 1., 0.], [-1., 0., 0.]]
        }
        ElementKind::Tri3 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]],
        _ => Vec::new(),
    }
}

/// Abaqus's shape functions of an element family at natural point `q`.
fn shape(kind: ElementKind, q: [f64; 3]) -> Vec<f64> {
    let nodes = natural(kind);
    let [x, y, z] = q;
    match kind {
        ElementKind::Hex8 => nodes.iter().map(|p| (1. + x * p[0]) * (1. + y * p[1]) * (1. + z * p[2]) / 8.).collect(),
        ElementKind::Hex20 => nodes
            .iter()
            .map(|p| {
                let (a, b, c) = (1. + x * p[0], 1. + y * p[1], 1. + z * p[2]);
                if p[0] == 0. {
                    (1. - x * x) * b * c / 4.
                } else if p[1] == 0. {
                    a * (1. - y * y) * c / 4.
                } else if p[2] == 0. {
                    a * b * (1. - z * z) / 4.
                } else {
                    a * b * c * (x * p[0] + y * p[1] + z * p[2] - 2.) / 8.
                }
            })
            .collect(),
        ElementKind::Tet4 => vec![1. - x - y - z, x, y, z],
        ElementKind::Tet10 => {
            let l = [1. - x - y - z, x, y, z];
            let mut n: Vec<f64> = l.iter().map(|li| li * (2. * li - 1.)).collect();
            for (a, b) in [(0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)] {
                n.push(4. * l[a] * l[b]);
            }
            n
        }
        ElementKind::Wedge6 => {
            let t = [1. - x - y, x, y];
            t.iter().map(|ti| ti * (1. - z) / 2.).chain(t.iter().map(|ti| ti * (1. + z) / 2.)).collect()
        }
        ElementKind::Quad4 => nodes.iter().map(|p| (1. + x * p[0]) * (1. + y * p[1]) / 4.).collect(),
        ElementKind::Quad8 => nodes
            .iter()
            .map(|p| {
                if p[0] == 0. {
                    (1. - x * x) * (1. + y * p[1]) / 2.
                } else if p[1] == 0. {
                    (1. + x * p[0]) * (1. - y * y) / 2.
                } else {
                    (1. + x * p[0]) * (1. + y * p[1]) * (x * p[0] + y * p[1] - 1.) / 4.
                }
            })
            .collect(),
        ElementKind::Tri3 => vec![1. - x - y, x, y],
        _ => Vec::new(),
    }
}

/// Columns ∂X/∂ξ, ∂X/∂η, ∂X/∂ζ of the mapping at `q` (central differences).
fn jacobian(kind: ElementKind, xyz: &[[f64; 3]], q: [f64; 3]) -> [[f64; 3]; 3] {
    let map = |q: [f64; 3]| {
        let n = shape(kind, q);
        let mut p = [0.0; 3];
        for (ni, x) in n.iter().zip(xyz) {
            for k in 0..3 {
                p[k] += ni * x[k];
            }
        }
        p
    };
    let h = 1e-6;
    [0, 1, 2].map(|d| {
        let (mut a, mut b) = (q, q);
        a[d] += h;
        b[d] -= h;
        let (pa, pb) = (map(a), map(b));
        [0, 1, 2].map(|k| (pa[k] - pb[k]) / (2. * h))
    })
}

fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
    [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]
}

fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
}

/// The smallest Jacobian determinant over the nodes and the centroid (solids), or the smallest
/// normal component along `up` (shells).
fn min_jacobian(kind: ElementKind, xyz: &[[f64; 3]], up: [f64; 3]) -> f64 {
    let nodes = natural(kind);
    let centroid = [0, 1, 2].map(|k| nodes.iter().map(|p| p[k]).sum::<f64>() / nodes.len() as f64);
    nodes
        .iter()
        .chain([&centroid])
        .map(|q| {
            let [a, b, c] = jacobian(kind, xyz, *q);
            if kind.dimension() == 3 { dot(cross(a, b), c) } else { dot(cross(a, b), up) }
        })
        .fold(f64::INFINITY, f64::min)
}

/// A one-element model of `kind` with a slanted, stretched, positive-volume placement of the
/// model's reference element (mid-side nodes at edge midpoints).
fn one_element(kind: ElementKind) -> Model {
    let corners: Vec<[f64; 3]> = match kind {
        ElementKind::Tri3 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]],
        ElementKind::Quad4 | ElementKind::Quad8 => vec![[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.]],
        ElementKind::Tet4 | ElementKind::Tet10 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]],
        ElementKind::Wedge6 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [0., 1., 1.]],
        _ => vec![[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]],
    };
    let mut xyz = corners.clone();
    xyz.resize(kind.node_count(), [0.0; 3]);
    for e in kind.edges() {
        if let [a, b, m] = **e {
            xyz[m] = [0, 1, 2].map(|i| (corners[a][i] + corners[b][i]) / 2.0);
        }
    }
    // A shear and stretch with a positive determinant, then a shift.
    let place = |p: [f64; 3]| [3. * p[0] + 0.5 * p[1] + 10., 2. * p[1] + 0.3 * p[2] - 4., 1.5 * p[2] + 0.2 * p[0] + 7.];
    let mut m = Model::default();
    for (i, p) in xyz.iter().enumerate() {
        m.set_node(i as u32 + 1, place(*p)).unwrap();
    }
    m.set_element(1, kind, (1..=kind.node_count() as u32).collect()).unwrap();
    m.add_set("all", elements([1])).unwrap();
    m.add_set("base", Set { nodes: vec![1], ..Set::default() }).unwrap();
    m.add_material("Steel", steel()).unwrap();
    let (section, formulation) = if kind.dimension() == 3 {
        (SectionKind::Solid {}, Formulation::default())
    } else {
        (SectionKind::Shell { thickness: 1.0, offset: 0.0, points: 5 }, Formulation { integration: Integration::Reduced, ..Formulation::default() })
    };
    m.add_section("s", Section { set: "all".into(), material: "Steel".into(), kind: section, formulation }).unwrap();
    m.add_boundary_condition("fix", BoundaryCondition::Fixed { target: Target::Set("base".into()) }).unwrap();
    m.add_step(step("Modes", Procedure::Frequency { modes: 1, min_hz: None, max_hz: None }, &["fix"], &[]), None).unwrap();
    m
}

#[test]
fn written_connectivity_has_a_positive_jacobian_in_abaqus_conventions() {
    let kinds = [ElementKind::Tet4, ElementKind::Tet10, ElementKind::Wedge6, ElementKind::Hex8, ElementKind::Hex20];
    let shells = [ElementKind::Tri3, ElementKind::Quad4, ElementKind::Quad8];
    for kind in kinds.into_iter().chain(shells) {
        let model = one_element(kind);
        let p = read(&abaqus::write(&model).unwrap());
        let (_, nodes) = &p.elements[&1];
        let xyz: Vec<[f64; 3]> = nodes.iter().map(|n| p.nodes[n]).collect();
        // A shell's positive normal is the side its face 1 (SPOS) looks to: the placed +Z here.
        let up = cross([3., 0., 0.2], [0.5, 2., 0.]);
        let j = min_jacobian(kind, &xyz, up);
        assert!(j > 0.0, "{kind:?}: smallest Jacobian {j}");

        // The test has teeth: a mirrored linear element shows up (and a misplaced quadratic
        // mid-side block, below).
        if kind.node_count() == natural(kind).len()
            && matches!(kind, ElementKind::Tet4 | ElementKind::Wedge6 | ElementKind::Hex8 | ElementKind::Tri3 | ElementKind::Quad4)
        {
            let mut mirrored = xyz.clone();
            match kind {
                ElementKind::Tet4 | ElementKind::Tri3 => mirrored.swap(1, 2),
                ElementKind::Wedge6 => {
                    mirrored.swap(1, 2);
                    mirrored.swap(4, 5);
                }
                ElementKind::Quad4 => mirrored.swap(1, 3),
                _ => {
                    mirrored.swap(1, 3);
                    mirrored.swap(5, 7);
                }
            }
            assert!(min_jacobian(kind, &mirrored, up) < 0.0, "{kind:?}: mirrored");
        }
    }
    // Hex20 in the order with the vertical edges' mid-sides before the top face's (Nastran's).
    let model = one_element(ElementKind::Hex20);
    let p = read(&abaqus::write(&model).unwrap());
    let mut xyz: Vec<[f64; 3]> = p.elements[&1].1.iter().map(|n| p.nodes[n]).collect();
    let (top, vertical) = (xyz[12..16].to_vec(), xyz[16..20].to_vec());
    xyz[12..16].copy_from_slice(&vertical);
    xyz[16..20].copy_from_slice(&top);
    assert!(min_jacobian(ElementKind::Hex20, &xyz, [0.; 3]) < 0.0, "a misplaced mid-side block is caught");
}

#[test]
fn model_faces_are_abaqus_faces() {
    // Each face's node positions (1-based, Abaqus element library, our transcription).
    let table: [(ElementKind, &[&[usize]]); 3] = [
        (ElementKind::Hex8, &[&[1, 2, 3, 4], &[5, 8, 7, 6], &[1, 5, 6, 2], &[2, 6, 7, 3], &[3, 7, 8, 4], &[4, 8, 5, 1]]),
        (ElementKind::Tet4, &[&[1, 2, 3], &[1, 4, 2], &[2, 4, 3], &[3, 4, 1]]),
        (ElementKind::Wedge6, &[&[1, 2, 3], &[4, 6, 5], &[1, 4, 5, 2], &[2, 5, 6, 3], &[3, 6, 4, 1]]),
    ];
    for (kind, faces) in table {
        assert_eq!(kind.faces().len(), faces.len());
        for (i, (ours, theirs)) in kind.faces().iter().zip(faces).enumerate() {
            let ours: BTreeSet<usize> = ours.iter().map(|p| p + 1).collect();
            let theirs: BTreeSet<usize> = theirs.iter().copied().collect();
            assert_eq!(ours, theirs, "{kind:?} face {}", i + 1);
            assert_eq!(abaqus::face_id(kind, i as u8 + 1), Some(format!("S{}", i + 1).as_str()));
        }
        assert_eq!(abaqus::face_id(kind, 0), None);
        assert_eq!(abaqus::face_id(kind, faces.len() as u8 + 1), None);
    }
    assert_eq!(abaqus::face_id(ElementKind::Quad4, 1), Some("SPOS"));
    assert_eq!(abaqus::face_id(ElementKind::Tri3, 2), Some("SNEG"));
    assert_eq!(abaqus::face_id(ElementKind::Beam2, 1), None);
}

// ---------------------------------------------------------------------------------------------
// Numbers and limits

#[test]
fn numbers_read_back_exactly() {
    for (v, s) in [(0.1, "0.1"), (210_000.0, "210000.0"), (7.85e-9, "7.85E-9"), (1e300, "1.0E300"), (-0.5, "-0.5"), (1e16, "1.0E16")] {
        assert_eq!(abaqus::number(v).unwrap(), s);
    }
    let mut x: u64 = 0x9E37_79B9_7F4A_7C15;
    for _ in 0..100_000 {
        x ^= x << 13;
        x ^= x >> 7;
        x ^= x << 17;
        let v = f64::from_bits(x);
        if !v.is_finite() {
            assert!(abaqus::number(v).is_err());
            continue;
        }
        let s = abaqus::number(v).unwrap();
        assert_eq!(s.parse::<f64>().unwrap().to_bits(), v.to_bits(), "{v:e} → {s}");
        assert!(s.len() <= 24, "{s}");
        assert!(!s.contains('e'), "{s}");
    }
    for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
        assert!(matches!(abaqus::number(bad), Err(Error::Unsupported { .. })));
    }
}

#[test]
fn largest_labels_longest_numbers_and_names_fit_a_line() {
    let mut m = Model::default();
    let base = modelcraft_model::MAX_LABEL - 30;
    let ugly = -1.234_567_890_123_456_7e-300;
    let mut xyz = one_element(ElementKind::Hex20).nodes().values().map(|n| n.xyz).collect::<Vec<_>>();
    xyz[0] = [ugly, ugly, ugly];
    for (i, p) in xyz.iter().enumerate() {
        m.set_node(base + i as u32, *p).unwrap();
    }
    m.set_element(modelcraft_model::MAX_LABEL, ElementKind::Hex20, (0..20).map(|i| base + i).collect()).unwrap();
    let long = format!("A{}", "x".repeat(modelcraft_model::MAX_NAME_LEN - 1));
    m.add_set(&long, Set { elements: vec![modelcraft_model::MAX_LABEL], nodes: (0..20).map(|i| base + i).collect(), ..Set::default() }).unwrap();
    m.add_material(&long, Material { elastic: Some(Elastic::Isotropic { young: ugly.abs(), poisson: 0.3 }), density: Some(1.0), fluid: None })
        .unwrap();
    m.add_section(&long, Section { set: long.clone(), material: long.clone(), kind: SectionKind::Solid {}, formulation: Formulation::default() })
        .unwrap();
    m.add_reference_point(&long, ReferencePoint { xyz: [ugly; 3], label: 0 }).unwrap();
    m.add_coupling(&long, Coupling { reference_point: long.clone(), set: long.clone(), kind: CouplingKind::Kinematic, dofs: Dofs::all() }).unwrap();
    m.add_boundary_condition(&long, BoundaryCondition::Fixed { target: Target::Set(long.clone()) }).unwrap();
    m.add_step(step(&long, Procedure::Frequency { modes: 100_000, min_hz: Some(ugly.abs()), max_hz: Some(1e300) }, &[&long], &[]), None).unwrap();
    let deck = abaqus::write(&m).unwrap();
    assert_lines_fit(&deck);
    let p = read(&deck);
    assert_eq!(p.elements[&modelcraft_model::MAX_LABEL].1.len(), 20, "a hex20 continues onto a second line");
    assert_eq!(p.nodes[&base], [ugly; 3]);
}

// ---------------------------------------------------------------------------------------------
// Step activity, load summing

fn block<'a>(deck: &'a str, step: &str) -> Vec<&'a str> {
    let start = format!("*STEP, NAME={step}");
    deck.lines().skip_while(|l| !l.starts_with(&start)).take_while(|l| *l != "*END STEP").filter(|l| !l.starts_with("**")).collect()
}

#[test]
fn steps_write_the_activity_table() {
    let deck = abaqus::write(&frequency()).unwrap();
    // `clamp` holds at zero in every step: the initial state.
    let initial: Vec<&str> = deck.lines().skip_while(|l| !l.starts_with("** Initial state")).skip(1).take(3).collect();
    assert_eq!(initial, ["*BOUNDARY", "** boundary condition clamp", "root, ENCASTRE"]);
    let twist = block(&deck, "Twist");
    assert!(twist.contains(&"*CLOAD, OP=NEW") && twist.contains(&"99999999, 4, 5000.0"), "{twist:?}");
    assert!(!twist.iter().any(|l| l.starts_with("*BOUNDARY")), "{twist:?}");
    // The frequency step adds `lift` to its base state and writes no loads.
    let modes = block(&deck, "Modes");
    assert!(modes.contains(&"*BOUNDARY") && modes.contains(&"99999999, 3, 3, 0.5"), "{modes:?}");
    assert!(!modes.iter().any(|l| l.contains("LOAD")), "{modes:?}");
    assert!(modes.contains(&"6, 1.0"), "{modes:?}");
    // `Lift` adds `lift` again: the frequency step's change did not carry on.
    let lift = block(&deck, "Lift");
    assert!(lift.contains(&"*BOUNDARY") && lift.contains(&"99999999, 3, 3, 0.5"), "{lift:?}");
    assert!(lift.contains(&"*DLOAD, OP=NEW") && lift.contains(&"all, BY, 0.001"), "{lift:?}");
    // `Release` drops `lift` and every load.
    let release = block(&deck, "Release");
    assert!(release.contains(&"*BOUNDARY, OP=NEW") && release.contains(&"root, ENCASTRE"), "{release:?}");
    for k in ["*CLOAD, OP=NEW", "*DLOAD, OP=NEW"] {
        let at = release.iter().position(|l| *l == k).unwrap_or_else(|| panic!("{k} in {release:?}"));
        assert!(release.get(at + 1).is_none_or(|l| l.starts_with('*')), "{k} restates nothing: {release:?}");
    }
    assert!(!release.contains(&"*DSLOAD, OP=NEW"), "no pressure before, none now");
    assert!(release.contains(&"U") && !release.contains(&"*ELEMENT OUTPUT"), "{release:?}");
}

#[test]
fn loads_sharing_nodes_elements_or_surfaces_add_up() {
    let mut m = plate();
    m.add_load("gust", Load::Pressure { set: "top".into(), magnitude: 0.02 }).unwrap();
    m.add_set("half", elements([1, 2])).unwrap();
    m.add_load("extra", Load::Gravity { acceleration: [0.0, 0.0, -9810.0], set: Some("half".into()) }).unwrap();
    m.add_load("p1", Load::Force { target: Target::Set("mid".into()), vector: [3.0, 0.0, 0.0], distribution: Distribution::PerNode }).unwrap();
    m.add_load("p2", Load::Force { target: Target::Set("mid".into()), vector: [1.0, 0.0, 0.0], distribution: Distribution::PerNode }).unwrap();
    for l in ["gust", "extra", "p1", "p2"] {
        m.set_load_active("Load", l, true).unwrap();
    }
    let deck = abaqus::write(&m).unwrap();
    let lines = block(&deck, "Load");
    assert!(lines.contains(&"top, P, 0.03"), "{lines:?}");
    assert!(lines.contains(&"5, 1, 4.0"), "{lines:?}");
    // Two gravity loads on elements 1 and 2: summed per element.
    assert!(lines.contains(&"1, GRAV, 19620.0, 0.0, 0.0, -1.0"), "{lines:?}");
    assert!(lines.contains(&"3, GRAV, 9810.0, 0.0, 0.0, -1.0"), "{lines:?}");
    assert!(!lines.iter().any(|l| l.starts_with("skin, GRAV")), "{lines:?}");
}

// ---------------------------------------------------------------------------------------------
// What the format cannot express, and models that are not ready

fn unsupported(m: &Model) -> String {
    match abaqus::write(m) {
        Err(Error::Unsupported { format: "Abaqus", what }) => what,
        other => panic!("expected unsupported, got {other:?}"),
    }
}

fn with_section(kind: ElementKind, section: SectionKind, integration: Integration) -> Model {
    let mut m = one_element(kind);
    let formulation = Formulation { integration, ..Formulation::default() };
    m.replace_section("s", Section { set: "all".into(), material: "Steel".into(), kind: section, formulation }).unwrap();
    m
}

#[test]
fn what_abaqus_cannot_express_is_named() {
    let shell = SectionKind::Shell { thickness: 1.0, offset: 0.0, points: 5 };
    let mut tri6 = Model::default();
    for (i, p) in [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0.5, 0., 0.], [0.5, 0.5, 0.], [0., 0.5, 0.]].iter().enumerate() {
        tri6.set_node(i as u32 + 1, *p).unwrap();
    }
    tri6.set_element(1, ElementKind::Tri6, (1..=6).collect()).unwrap();
    tri6.add_set("all", elements([1])).unwrap();
    tri6.add_material("Steel", steel()).unwrap();
    tri6.add_section("s", Section { set: "all".into(), material: "Steel".into(), kind: shell, formulation: Formulation::default() }).unwrap();
    tri6.add_step(step("Modes", Procedure::Frequency { modes: 1, min_hz: None, max_hz: None }, &[], &[]), None).unwrap();
    assert!(unsupported(&tri6).contains("tri6 shells"));

    assert!(unsupported(&with_section(ElementKind::Quad8, shell, Integration::Full)).contains("full integration of quad8"));
    assert!(unsupported(&with_section(ElementKind::Tet10, SectionKind::Solid {}, Integration::Reduced)).contains("reduced integration of tet10"));
    assert!(unsupported(&with_section(ElementKind::Wedge6, SectionKind::Solid {}, Integration::Reduced)).contains("wedge6"));
    // A linear tetrahedron integrates at one point either way.
    assert!(abaqus::write(&with_section(ElementKind::Tet4, SectionKind::Solid {}, Integration::Reduced)).unwrap().contains("TYPE=C3D4\n"));

    let mut case = cantilever();
    case.add_set("ROOT", Set { nodes: vec![1], ..Set::default() }).unwrap();
    let what = unsupported(&case);
    assert!(what.contains("sets `ROOT` and `root`") || what.contains("sets `root` and `ROOT`"), "{what}");

    let mut ortho = plate();
    let carbon = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] };
    ortho.replace_material("Aluminium", Material { elastic: Some(carbon), density: Some(1.6e-9), fluid: None }).unwrap();
    assert!(unsupported(&ortho).contains("orthotropic material `Aluminium` in a shell section"));

    let mut loose = frequency();
    let dofs = Dofs::try_from(vec![1, 2, 4, 5, 6]).unwrap();
    loose.replace_coupling("tie", Coupling { reference_point: "Tip".into(), set: "tip".into(), kind: CouplingKind::Distributing, dofs }).unwrap();
    assert!(unsupported(&loose).contains("coupling `tie`"));

    // Several problems are named together.
    let mut both = with_section(ElementKind::Quad8, shell, Integration::Full);
    both.add_set("ALL", Set { nodes: vec![1], ..Set::default() }).unwrap();
    assert!(abaqus::write(&both).unwrap_err().to_string().contains("full integration"));
}

#[test]
fn orthotropic_solids_beams_and_rods_are_written() {
    let mut m = cantilever();
    let carbon = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] };
    m.replace_material("Steel", Material { elastic: Some(carbon), density: None, fluid: None }).unwrap();
    let deck = abaqus::write(&m).unwrap();
    assert!(deck.contains("*ELASTIC, TYPE=ENGINEERING CONSTANTS\n140000.0, 10000.0, 10000.0, 0.3, 0.3, 0.4, 5000.0, 5000.0\n3800.0\n"), "{deck}");

    // A frame of one beam and one rod.
    let mut f = Model::default();
    for (i, p) in [[0., 0., 0.], [1000., 0., 0.], [1000., 0., 500.]].iter().enumerate() {
        f.set_node(i as u32 + 1, *p).unwrap();
    }
    f.set_element(1, ElementKind::Beam2, vec![1, 2]).unwrap();
    f.set_element(2, ElementKind::Rod2, vec![2, 3]).unwrap();
    f.add_set("beam", elements([1])).unwrap();
    f.add_set("rod", elements([2])).unwrap();
    f.add_set("base", Set { nodes: vec![1, 3], ..Set::default() }).unwrap();
    f.add_material("Steel", steel()).unwrap();
    let i = BeamProfile::I { height: 100.0, width: 50.0, web: 5.0, flange: 8.0 };
    let beam = SectionKind::Beam { profile: i, orientation: [0.0, 1.0, 0.0] };
    f.add_section("b", Section { set: "beam".into(), material: "Steel".into(), kind: beam, formulation: Formulation::default() }).unwrap();
    let rod = SectionKind::Rod { area: 50.0 };
    f.add_section("r", Section { set: "rod".into(), material: "Steel".into(), kind: rod, formulation: Formulation::default() }).unwrap();
    f.add_boundary_condition("pin", BoundaryCondition::Fixed { target: Target::Set("base".into()) }).unwrap();
    f.add_step(step("Modes", Procedure::Frequency { modes: 3, min_hz: None, max_hz: Some(500.0) }, &["pin"], &[]), None).unwrap();
    let deck = abaqus::write(&f).unwrap();
    for expected in [
        "*ELEMENT, TYPE=B31\n1, 1, 2\n",
        "*ELEMENT, TYPE=T3D2\n2, 2, 3\n",
        "*BEAM SECTION, SECTION=I, ELSET=beam, MATERIAL=Steel\n50.0, 100.0, 50.0, 50.0, 8.0, 8.0, 5.0\n0.0, 1.0, 0.0\n",
        "*SOLID SECTION, ELSET=rod, MATERIAL=Steel\n50.0\n",
        "3, , 500.0\n",
    ] {
        assert!(deck.contains(expected), "{expected}\n{deck}");
    }
}

#[test]
fn a_model_that_is_not_ready_is_refused() {
    let mut m = cantilever();
    m.remove_section("beam").unwrap();
    assert!(matches!(abaqus::write(&m), Err(Error::Model(_))));
    assert!(matches!(abaqus::write(&Model::default()), Err(Error::Model(_))));
}