forked from 1b-agent/fieldcraft-upstream
jyvikoandClaude Opus 5.5
3a5a8213a1
A five-node pyramid, the conforming transition between a hexahedron's face and tetrahedra: Bedrosian's rational shape functions with a collapsed (Duffy) Gauss rule, a certified Jacobian on the base square with exact volume, a scaled Jacobian, C3D5 / CPYRAM / LS-DYNA degenerate solid, VTK_PYRAMID and the binary mesh code 13. A distorted hex8-pyramid5-tet4 transition mesh passes the patch test and uniaxial tension exactly. Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo
871 lines
46 KiB
Rust
871 lines
46 KiB
Rust
//! 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 fieldcraft_decks::{Error, abaqus};
|
||
use fieldcraft_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), ..Material::default() }
|
||
}
|
||
|
||
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), ..Material::default() };
|
||
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("FIELDCRAFT_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 FIELDCRAFT_BLESS=1 to write it)", path.display()));
|
||
assert!(expected == deck, "{name}.inp differs from the written deck (FIELDCRAFT_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([¢roid])
|
||
.map(|q| {
|
||
let [a, b, c] = jacobian(kind, xyz, *q);
|
||
if kind.dimension() == 3 { dot(cross(a, b), c) } else { dot(cross(a, b), up) }
|
||
})
|
||
.fold(f64::INFINITY, f64::min)
|
||
}
|
||
|
||
/// A one-element model of `kind` with a slanted, stretched, positive-volume placement of the
|
||
/// model's reference element (mid-side nodes at edge midpoints).
|
||
fn one_element(kind: ElementKind) -> Model {
|
||
let corners: Vec<[f64; 3]> = match kind {
|
||
ElementKind::Tri3 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]],
|
||
ElementKind::Quad4 | ElementKind::Quad8 => vec![[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.]],
|
||
ElementKind::Tet4 | ElementKind::Tet10 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]],
|
||
ElementKind::Wedge6 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [0., 1., 1.]],
|
||
_ => vec![[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]],
|
||
};
|
||
let mut xyz = corners.clone();
|
||
xyz.resize(kind.node_count(), [0.0; 3]);
|
||
for e in kind.edges() {
|
||
if let [a, b, m] = **e {
|
||
xyz[m] = [0, 1, 2].map(|i| (corners[a][i] + corners[b][i]) / 2.0);
|
||
}
|
||
}
|
||
// A shear and stretch with a positive determinant, then a shift.
|
||
let place = |p: [f64; 3]| [3. * p[0] + 0.5 * p[1] + 10., 2. * p[1] + 0.3 * p[2] - 4., 1.5 * p[2] + 0.2 * p[0] + 7.];
|
||
let mut m = Model::default();
|
||
for (i, p) in xyz.iter().enumerate() {
|
||
m.set_node(i as u32 + 1, place(*p)).unwrap();
|
||
}
|
||
m.set_element(1, kind, (1..=kind.node_count() as u32).collect()).unwrap();
|
||
m.add_set("all", elements([1])).unwrap();
|
||
m.add_set("base", Set { nodes: vec![1], ..Set::default() }).unwrap();
|
||
m.add_material("Steel", steel()).unwrap();
|
||
let (section, formulation) = if kind.dimension() == 3 {
|
||
(SectionKind::Solid {}, Formulation::default())
|
||
} else {
|
||
(SectionKind::Shell { thickness: 1.0, offset: 0.0, points: 5 }, Formulation { integration: Integration::Reduced, ..Formulation::default() })
|
||
};
|
||
m.add_section("s", Section { set: "all".into(), material: "Steel".into(), kind: section, formulation }).unwrap();
|
||
m.add_boundary_condition("fix", BoundaryCondition::Fixed { target: Target::Set("base".into()) }).unwrap();
|
||
m.add_step(step("Modes", Procedure::Frequency { modes: 1, min_hz: None, max_hz: None }, &["fix"], &[]), None).unwrap();
|
||
m
|
||
}
|
||
|
||
#[test]
|
||
fn written_connectivity_has_a_positive_jacobian_in_abaqus_conventions() {
|
||
let kinds = [ElementKind::Tet4, ElementKind::Tet10, ElementKind::Wedge6, ElementKind::Hex8, ElementKind::Hex20];
|
||
let shells = [ElementKind::Tri3, ElementKind::Quad4, ElementKind::Quad8];
|
||
for kind in kinds.into_iter().chain(shells) {
|
||
let model = one_element(kind);
|
||
let p = read(&abaqus::write(&model).unwrap());
|
||
let (_, nodes) = &p.elements[&1];
|
||
let xyz: Vec<[f64; 3]> = nodes.iter().map(|n| p.nodes[n]).collect();
|
||
// A shell's positive normal is the side its face 1 (SPOS) looks to: the placed +Z here.
|
||
let up = cross([3., 0., 0.2], [0.5, 2., 0.]);
|
||
let j = min_jacobian(kind, &xyz, up);
|
||
assert!(j > 0.0, "{kind:?}: smallest Jacobian {j}");
|
||
|
||
// The test has teeth: a mirrored linear element shows up (and a misplaced quadratic
|
||
// mid-side block, below).
|
||
if kind.node_count() == natural(kind).len()
|
||
&& matches!(kind, ElementKind::Tet4 | ElementKind::Wedge6 | ElementKind::Hex8 | ElementKind::Tri3 | ElementKind::Quad4)
|
||
{
|
||
let mut mirrored = xyz.clone();
|
||
match kind {
|
||
ElementKind::Tet4 | ElementKind::Tri3 => mirrored.swap(1, 2),
|
||
ElementKind::Wedge6 => {
|
||
mirrored.swap(1, 2);
|
||
mirrored.swap(4, 5);
|
||
}
|
||
ElementKind::Quad4 => mirrored.swap(1, 3),
|
||
_ => {
|
||
mirrored.swap(1, 3);
|
||
mirrored.swap(5, 7);
|
||
}
|
||
}
|
||
assert!(min_jacobian(kind, &mirrored, up) < 0.0, "{kind:?}: mirrored");
|
||
}
|
||
}
|
||
// Hex20 in the order with the vertical edges' mid-sides before the top face's (Nastran's).
|
||
let model = one_element(ElementKind::Hex20);
|
||
let p = read(&abaqus::write(&model).unwrap());
|
||
let mut xyz: Vec<[f64; 3]> = p.elements[&1].1.iter().map(|n| p.nodes[n]).collect();
|
||
let (top, vertical) = (xyz[12..16].to_vec(), xyz[16..20].to_vec());
|
||
xyz[12..16].copy_from_slice(&vertical);
|
||
xyz[16..20].copy_from_slice(&top);
|
||
assert!(min_jacobian(ElementKind::Hex20, &xyz, [0.; 3]) < 0.0, "a misplaced mid-side block is caught");
|
||
}
|
||
|
||
#[test]
|
||
fn model_faces_are_abaqus_faces() {
|
||
// Each face's node positions (1-based, Abaqus element library, our transcription).
|
||
let table: [(ElementKind, &[&[usize]]); 4] = [
|
||
(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]]),
|
||
(ElementKind::Pyramid5, &[&[1, 2, 3, 4], &[1, 5, 2], &[2, 5, 3], &[3, 5, 4], &[4, 5, 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 = fieldcraft_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(fieldcraft_model::MAX_LABEL, ElementKind::Hex20, (0..20).map(|i| base + i).collect()).unwrap();
|
||
let long = format!("A{}", "x".repeat(fieldcraft_model::MAX_NAME_LEN - 1));
|
||
m.add_set(&long, Set { elements: vec![fieldcraft_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), ..Material::default() },
|
||
)
|
||
.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[&fieldcraft_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), ..Material::default() }).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, ..Material::default() }).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(_))));
|
||
}
|
||
|
||
/// The cantilever with its steel given each plasticity in turn.
|
||
fn plastic_cantilever(p: fieldcraft_model::Plastic) -> Model {
|
||
let mut m = cantilever();
|
||
m.replace_material("Steel", Material { plastic: Some(p), ..steel() }).unwrap();
|
||
m
|
||
}
|
||
|
||
/// Plasticity is written as `*PLASTIC` (and Johnson–Cook's rate and heat keywords); damage and
|
||
/// erosion, failure models of Abaqus/Explicit, are reported.
|
||
#[test]
|
||
fn plasticity_is_written_and_failure_reported() {
|
||
use fieldcraft_model::{Damage, Erosion, Hardening, HardeningPoint, Plastic};
|
||
let section = |deck: &str| deck.split("*MATERIAL, NAME=Steel").nth(1).and_then(|s| s.split("*SOLID SECTION").next()).unwrap().to_string();
|
||
let deck =
|
||
abaqus::write(&plastic_cantilever(Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 1000.0, hardening: Hardening::Isotropic }))
|
||
.unwrap();
|
||
assert!(section(&deck).contains("*PLASTIC\n250.0, 0.0\n10250.0, 10.0\n"), "{}", section(&deck));
|
||
let deck =
|
||
abaqus::write(&plastic_cantilever(Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 1000.0, hardening: Hardening::Kinematic }))
|
||
.unwrap();
|
||
assert!(section(&deck).contains("*PLASTIC, HARDENING=KINEMATIC\n250.0, 0.0\n1250.0, 1.0\n"), "{}", section(&deck));
|
||
let table = Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: 200.0 }, HardeningPoint { strain: 0.1, stress: 300.0 }] };
|
||
let deck = abaqus::write(&plastic_cantilever(table)).unwrap();
|
||
assert!(section(&deck).contains("*PLASTIC\n200.0, 0.0\n300.0, 0.1\n"), "{}", section(&deck));
|
||
let jc = Plastic::JohnsonCook {
|
||
a: 90.0,
|
||
b: 292.0,
|
||
n: 0.31,
|
||
c: 0.025,
|
||
m: 1.09,
|
||
reference_strain_rate: 1.0,
|
||
room_temperature: 294.0,
|
||
melt_temperature: 1356.0,
|
||
specific_heat: Some(3.83e8),
|
||
taylor_quinney: 0.9,
|
||
};
|
||
let deck = abaqus::write(&plastic_cantilever(jc.clone())).unwrap();
|
||
let s = section(&deck);
|
||
assert!(
|
||
s.contains("*PLASTIC, HARDENING=JOHNSON COOK\n90.0, 292.0, 0.31, 1.09, 1356.0, 294.0\n*RATE DEPENDENT, TYPE=JOHNSON COOK\n0.025, 1.0\n*SPECIFIC HEAT\n383000000.0\n*INELASTIC HEAT FRACTION\n0.9\n"),
|
||
"{s}"
|
||
);
|
||
let mut failing = plastic_cantilever(jc);
|
||
let mut steel = failing.material("Steel").unwrap().clone();
|
||
steel.damage = Some(Damage::JohnsonCook { d1: 0.54, d2: 4.89, d3: -3.03, d4: 0.014, d5: 1.12 });
|
||
steel.erosion = Some(Erosion { plastic_strain: 1.0 });
|
||
failing.replace_material("Steel", steel).unwrap();
|
||
let Err(Error::Unsupported { what, .. }) = abaqus::write(&failing) else { panic!("damage is not written") };
|
||
assert!(what.contains("material `Steel`: damage and erosion"), "{what}");
|
||
}
|
||
|
||
/// A nonlinear static step is a general step with NLGEOM, its increments and its line search; a
|
||
/// neo-Hookean material is then `*HYPERELASTIC` and a linear static step beside it is refused.
|
||
#[test]
|
||
fn nonlinear_steps_and_hyperelastic_materials_are_written() {
|
||
use fieldcraft_model::{Hyperelastic, NonlinearStatic};
|
||
let mut m = cantilever();
|
||
let controls = NonlinearStatic {
|
||
initial_increment: 0.05,
|
||
min_increment: 1e-4,
|
||
max_increment: 0.5,
|
||
max_increments: 40,
|
||
line_search: true,
|
||
..NonlinearStatic::default()
|
||
};
|
||
m.replace_step(step("Bend", Procedure::NonlinearStatic(controls.clone()), &["clamp"], &["end"])).unwrap();
|
||
let deck = abaqus::write(&m).unwrap();
|
||
let bend = block(&deck, "Bend");
|
||
assert_eq!(bend.first(), Some(&"*STEP, NAME=Bend, NLGEOM=YES, INC=40"), "{bend:?}");
|
||
let at = bend.iter().position(|l| *l == "*STATIC").unwrap();
|
||
assert_eq!(bend.get(at + 1..at + 4), Some(&["0.05, 1.0, 0.0001, 0.5", "*CONTROLS, PARAMETERS=LINE SEARCH", "5,"][..]), "{bend:?}");
|
||
assert!(deck.contains("*ELASTIC, TYPE=ISOTROPIC\n210000.0, 0.3\n"), "St Venant–Kirchhoff is the elastic material");
|
||
|
||
let small = NonlinearStatic { nonlinear_geometry: false, line_search: false, ..controls };
|
||
m.replace_step(step("Bend", Procedure::NonlinearStatic(small), &["clamp"], &["end"])).unwrap();
|
||
m.replace_material("Steel", Material { hyperelastic: Some(Hyperelastic::NeoHookean {}), ..steel() }).unwrap();
|
||
let deck = abaqus::write(&m).unwrap();
|
||
let bend = block(&deck, "Bend");
|
||
assert_eq!(bend.first(), Some(&"*STEP, NAME=Bend, NLGEOM=NO, INC=40"), "{bend:?}");
|
||
assert!(!bend.iter().any(|l| l.starts_with("*CONTROLS")), "{bend:?}");
|
||
// μ = E/(2(1 + ν)) = 80769.23…, κ = E/(3(1 − 2ν)) = 175000: C10 = μ/2, D1 = 2/κ.
|
||
let (c10, d1) = (210_000.0 / 2.6 / 2.0, 2.0 / 175_000.0);
|
||
let line = deck.split("*HYPERELASTIC, NEO HOOKE\n").nth(1).and_then(|s| s.lines().next()).unwrap();
|
||
let values: Vec<f64> = line.split(", ").map(|v| v.parse().unwrap()).collect();
|
||
assert!((values[0] - c10).abs() < 1e-9 * c10 && (values[1] - d1).abs() < 1e-9 * d1, "{line}");
|
||
assert!(!deck.contains("*ELASTIC"), "{deck}");
|
||
|
||
// Mooney–Rivlin and Yeoh: C10 = μ/2 − C01, C01, D1 and C10, C20, C30, D1.
|
||
m.replace_material("Steel", Material { hyperelastic: Some(Hyperelastic::MooneyRivlin { c01: 10000.0 }), ..steel() }).unwrap();
|
||
let deck = abaqus::write(&m).unwrap();
|
||
let line = deck.split("*HYPERELASTIC, MOONEY-RIVLIN\n").nth(1).and_then(|s| s.lines().next()).unwrap();
|
||
let values: Vec<f64> = line.split(", ").map(|v| v.parse().unwrap()).collect();
|
||
assert!((values[0] - (c10 - 10000.0)).abs() < 1e-9 * c10 && values[1] == 10000.0 && (values[2] - d1).abs() < 1e-9 * d1, "{line}");
|
||
m.replace_material("Steel", Material { hyperelastic: Some(Hyperelastic::Yeoh { c20: -50.0, c30: 2.0 }), ..steel() }).unwrap();
|
||
let deck = abaqus::write(&m).unwrap();
|
||
let line = deck.split("*HYPERELASTIC, YEOH\n").nth(1).and_then(|s| s.lines().next()).unwrap();
|
||
let values: Vec<f64> = line.split(", ").map(|v| v.parse().unwrap()).collect();
|
||
assert!((values[0] - c10).abs() < 1e-9 * c10 && values[1..3] == [-50.0, 2.0] && (values[3] - d1).abs() < 1e-9 * d1, "{line}");
|
||
m.replace_material("Steel", Material { hyperelastic: Some(Hyperelastic::NeoHookean {}), ..steel() }).unwrap();
|
||
|
||
m.add_step(step("Linear", Procedure::Static {}, &["clamp"], &["end"]), None).unwrap();
|
||
assert!(unsupported(&m).contains("step `Linear`: a linear static step beside nonlinear steps of hyperelastic material `Steel`"));
|
||
m.remove_step("Bend").unwrap();
|
||
assert!(abaqus::write(&m).unwrap().contains("*ELASTIC, TYPE=ISOTROPIC\n"), "with linear steps only, the elastic constants are the whole of it");
|
||
}
|
||
|
||
/// A plastic material in nonlinear steps is `*PLASTIC` in a general step; what Abaqus would
|
||
/// read differently is refused: a load path, and several general steps (Abaqus carries the
|
||
/// plastic state between them).
|
||
#[test]
|
||
fn nonlinear_plastic_steps_and_what_is_refused() {
|
||
use fieldcraft_model::{Hardening, NonlinearStatic, Plastic};
|
||
let mut m = cantilever();
|
||
let plastic = Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 2000.0, hardening: Hardening::Isotropic };
|
||
m.replace_material("Steel", Material { plastic: Some(plastic), ..steel() }).unwrap();
|
||
m.replace_step(step("Bend", Procedure::NonlinearStatic(NonlinearStatic::default()), &["clamp"], &["end"])).unwrap();
|
||
let deck = abaqus::write(&m).unwrap();
|
||
assert!(deck.contains("*PLASTIC\n250.0, 0.0\n") && deck.contains("NLGEOM=YES"), "{deck}");
|
||
let path = NonlinearStatic { load_path: Some(vec![[0.0, 0.0], [0.5, 1.0], [1.0, 0.0]]), ..NonlinearStatic::default() };
|
||
m.replace_step(step("Bend", Procedure::NonlinearStatic(path), &["clamp"], &["end"])).unwrap();
|
||
assert!(unsupported(&m).contains("step `Bend`: a load path"));
|
||
m.replace_step(step("Bend", Procedure::NonlinearStatic(NonlinearStatic::default()), &["clamp"], &["end"])).unwrap();
|
||
m.add_step(step("Again", Procedure::NonlinearStatic(NonlinearStatic::default()), &["clamp"], &["end"]), None).unwrap();
|
||
assert!(unsupported(&m).contains("plastic material `Steel` in a deck of 2 general steps"));
|
||
}
|
||
|
||
/// A buckling step is `*BUCKLE` with its loads as the perturbation loads; after a general step
|
||
/// (Abaqus's base state would be that step's) and with a shift it is refused by name.
|
||
#[test]
|
||
fn buckling_steps() {
|
||
let mut m = cantilever_mesh();
|
||
m.add_load("end", Load::Force { target: Target::Set("tip".into()), vector: [-100.0, 0.0, 0.0], distribution: Distribution::Total }).unwrap();
|
||
m.add_step(step("Buckle", Procedure::Buckle { modes: 4, shift: None }, &["clamp"], &["end"]), None).unwrap();
|
||
let deck = abaqus::write(&m).unwrap();
|
||
assert!(deck.contains("*BUCKLE, EIGENSOLVER=LANCZOS\n4\n") && deck.contains("*CLOAD"), "{deck}");
|
||
let mut shifted = m.clone();
|
||
shifted.replace_step(step("Buckle", Procedure::Buckle { modes: 4, shift: Some(2.0) }, &["clamp"], &["end"])).unwrap();
|
||
assert!(abaqus::write(&shifted).unwrap_err().to_string().contains("shift"));
|
||
m.add_step(step("Bend", Procedure::Static {}, &["clamp"], &["end"]), Some("Buckle")).unwrap();
|
||
assert!(abaqus::write(&m).unwrap_err().to_string().contains("after a general step"));
|
||
assert!(fieldcraft_decks::lsdyna::write(&m).is_err() && fieldcraft_decks::nastran::write(&m).is_err());
|
||
}
|
||
|
||
/// An arc-length step is `*STATIC, RIKS` with the arcs, the maximum load factor and the
|
||
/// displacement that end it; a spherical constraint is refused, and LS-DYNA refuses the step.
|
||
#[test]
|
||
fn arc_length_steps() {
|
||
use fieldcraft_model::{ArcConstraint, ArcLength, DisplacementLimit, NonlinearStatic};
|
||
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();
|
||
let arc = ArcLength {
|
||
max_load_factor: 3.0,
|
||
displacement_limit: Some(DisplacementLimit { node: 2, direction: 3, displacement: -5.0 }),
|
||
..ArcLength::default()
|
||
};
|
||
let settings = NonlinearStatic { arc_length: Some(arc.clone()), ..NonlinearStatic::default() };
|
||
m.add_step(step("Snap", Procedure::NonlinearStatic(settings), &["clamp"], &["end"]), None).unwrap();
|
||
let deck = abaqus::write(&m).unwrap();
|
||
assert!(deck.contains("*STATIC, RIKS\n0.1, 1.0, 1.0E-5, 1.0, 3.0, 2, 3, -5.0\n"), "{deck}");
|
||
assert!(fieldcraft_decks::lsdyna::write(&m).is_err());
|
||
let spherical = NonlinearStatic { arc_length: Some(ArcLength { constraint: ArcConstraint::Spherical, ..arc }), ..NonlinearStatic::default() };
|
||
m.replace_step(step("Snap", Procedure::NonlinearStatic(spherical), &["clamp"], &["end"])).unwrap();
|
||
assert!(abaqus::write(&m).is_err());
|
||
}
|