diff --git a/ROADMAP.md b/ROADMAP.md index 0dfdbb9..0499320 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -21,7 +21,7 @@ hex meshes today and it gives the best elements where stresses are read. | M0 Foundation | Workspace, quality gates, FE model, command engine with rollback and panic guard, CLI, geometry boundary on SolidCraft | done | | M1 Geometry in and surface preparation | Open `.solidcraft` designs with persistent face names, assemblies placed (done); face-tagged 3MF (done); test corpus, `cargo xtask corpus` (done); exact predicates (done); quality metrics and certified validity gate (done); weld, closed/manifold/orientation and self-intersection checks, feature edges (`modelcraft-mesh::surface`, `mesh.surface.check`, done) | in progress | | M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity; the `modelcraft_mesh::size::SizeField` trait and edge-division helpers are in); Bowyer–Watson Delaunay (exact, symbolic perturbation, verified; `modelcraft-mesh::delaunay`, done); robust CDT boundary recovery (Diazzi et al. 2023); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress | -| M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; writers for Abaqus input, Nastran SOL 101/103 and LS-DYNA implicit; the `.modelcraft` file | planned | +| M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input, Nastran SOL 101/103 and LS-DYNA implicit; the `.modelcraft` file | in progress | | M4 Own solver and first results | Linear static and modal solver (tet4/10, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions; deformed shape, contours, mode shapes, probes, reaction sums | planned | | M5 Structured hex | 2D quad meshing of faces (`modelcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | | M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned | diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs new file mode 100644 index 0000000..a6b8880 --- /dev/null +++ b/crates/engine/src/cmd/define.rs @@ -0,0 +1,385 @@ +//! Model definition commands: create, edit and delete sets, materials, sections, reference +//! points, couplings, boundary conditions, loads and steps, and the step activity table. +//! +//! Each definition's parameters are its JSON form in the model (the same form a `.modelcraft` +//! file stores), plus its `name`. `create` fails when the name is taken, `edit` replaces an +//! existing definition, `delete` refuses while another definition refers to it. + +use modelcraft_model::{BoundaryCondition, Coupling, Load, Material, Model, ReferencePoint, Section, Set, Step}; +use serde::de::DeserializeOwned; +use serde_json::{Value, json}; + +use super::{name, opt, req}; +use crate::{CommandSpec, EngineError, Result, Session}; + +/// The parameters other than `name`, read as a definition's JSON form. +fn fields(p: &Value, what: &str) -> Result { + let mut obj = p.as_object().cloned().ok_or_else(|| EngineError::Params(format!("the parameters of a {what} are a JSON object")))?; + obj.remove("name"); + serde_json::from_value(Value::Object(obj)).map_err(|e| EngineError::Params(format!("{what}: {e}"))) +} + +/// Create, edit and delete commands for one kind of definition. +macro_rules! definition_commands { + ($what:literal, $ty:ty, $id:literal, $label:literal, $doc:literal, $add:ident, $replace:ident, $remove:ident, + $CREATE:ident, $EDIT:ident, $DELETE:ident) => { + pub const $CREATE: CommandSpec = CommandSpec { + id: concat!("model.", $id, ".create"), + label: concat!("New ", $label), + params: concat!($doc, " — fails if the name is taken"), + run: |s: &mut Session, p: &Value| { + let n = name(p)?; + s.model.$add(n, fields::<$ty>(p, $what)?)?; + Ok(json!({"name": n})) + }, + }; + pub const $EDIT: CommandSpec = CommandSpec { + id: concat!("model.", $id, ".edit"), + label: concat!("Edit ", $label), + params: concat!($doc, " — replaces the existing definition of that name"), + run: |s: &mut Session, p: &Value| { + let n = name(p)?; + s.model.$replace(n, fields::<$ty>(p, $what)?)?; + Ok(json!({"name": n})) + }, + }; + pub const $DELETE: CommandSpec = CommandSpec { + id: concat!("model.", $id, ".delete"), + label: concat!("Delete ", $label), + params: r#"{"name": "…"} — refused while another definition refers to it"#, + run: |s: &mut Session, p: &Value| { + s.model.$remove(name(p)?)?; + Ok(Value::Null) + }, + }; + }; +} + +definition_commands!( + "set", + Set, + "set", + "Set", + r#"{"name": "root", "geometry": [{"body": "Body1", "face": "F2:top"}, {"body": "Body1", "edge": "…"}, {"body": "Body1", "vertex": "…"}, {"body": "Body1"}], "nodes": [labels], "elements": [labels], "faces": [{"element": 1, "face": 3}]} — any of the member lists; geometry members resolve against the mesh's geometry tags after every remesh (a face name also selects its split pieces `name#0`, `name#1`, …)"#, + add_set, + replace_set, + remove_set, + SET_CREATE, + SET_EDIT, + SET_DELETE +); + +definition_commands!( + "material", + Material, + "material", + "Material", + r#"{"name": "Steel", "elastic": {"type": "isotropic", "young": 210000, "poisson": 0.3} | {"type": "engineering_constants", "young": [E1, E2, E3], "poisson": [ν12, ν13, ν23], "shear": [G12, G13, G23]}, "density": 7.85e-9} — in the model's unit system (model.units), never converted"#, + add_material, + replace_material, + remove_material, + MATERIAL_CREATE, + MATERIAL_EDIT, + MATERIAL_DELETE +); + +definition_commands!( + "section", + Section, + "section", + "Section", + r#"{"name": "web", "set": "all", "material": "Steel", "kind": {"type": "solid"} | {"type": "shell", "thickness": 2, "offset": 0 (-0.5..=0.5 of the thickness), "points": 5} | {"type": "beam", "profile": {"type": "rectangle", "width", "height"} | {"type": "circle", "radius"} | {"type": "pipe", "radius", "thickness"} | {"type": "i", "height", "width", "web", "flange"} | {"type": "box", "width", "height", "thickness"}, "orientation": [x, y, z]} | {"type": "rod", "area": 10}, "formulation": {"integration": "full|reduced", "incompatible_modes": false, "hybrid": false}}"#, + add_section, + replace_section, + remove_section, + SECTION_CREATE, + SECTION_EDIT, + SECTION_DELETE +); + +definition_commands!( + "reference point", + ReferencePoint, + "reference_point", + "Reference Point", + r#"{"name": "RP", "xyz": [x, y, z], "label": 0} — label: its node label in decks (default: the highest free label, kept through remeshes)"#, + add_reference_point, + replace_reference_point, + remove_reference_point, + REFERENCE_POINT_CREATE, + REFERENCE_POINT_EDIT, + REFERENCE_POINT_DELETE +); + +definition_commands!( + "coupling", + Coupling, + "coupling", + "Coupling", + r#"{"name": "pad", "reference_point": "RP", "set": "hole", "kind": "kinematic|distributing", "dofs": [1, 2, 3, 4, 5, 6] (default all)}"#, + add_coupling, + replace_coupling, + remove_coupling, + COUPLING_CREATE, + COUPLING_EDIT, + COUPLING_DELETE +); + +definition_commands!( + "boundary condition", + BoundaryCondition, + "boundary_condition", + "Boundary Condition", + r#"{"name": "clamp", "type": "fixed", "target": {"set": "root"} | {"reference_point": "RP"}} | {"name", "type": "displacement", "target", "values": [u1, u2, u3, ur1, ur2, ur3] (null leaves a degree of freedom free)} | {"name", "type": "symmetry", "target", "plane": "x|y|z" (the plane's normal)} — active in the steps that list it (model.step.activate)"#, + add_boundary_condition, + replace_boundary_condition, + remove_boundary_condition, + BOUNDARY_CONDITION_CREATE, + BOUNDARY_CONDITION_EDIT, + BOUNDARY_CONDITION_DELETE +); + +definition_commands!( + "load", + Load, + "load", + "Load", + r#"{"name": "push", "type": "force|moment", "target": {"set": "…"} | {"reference_point": "…"}, "vector": [x, y, z], "distribution": "total" (default: the vector is the total, split by area over the set's faces, or equally over its nodes) | "per_node" (at every node: depends on the mesh)} | {"name", "type": "traction", "set", "vector" (per area, on the set's faces)} | {"name", "type": "pressure", "set", "magnitude" (positive pushes in)} | {"name", "type": "gravity", "acceleration": [x, y, z], "set" (optional, default the whole model)} | {"name", "type": "body_force", "set", "vector" (per volume)} — active in the steps that list it (model.step.activate)"#, + add_load, + replace_load, + remove_load, + LOAD_CREATE, + LOAD_EDIT, + LOAD_DELETE +); + +pub const STEP_CREATE: CommandSpec = CommandSpec { + id: "model.step.create", + label: "New Step", + params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, + run: step_create, +}; + +pub const STEP_EDIT: CommandSpec = CommandSpec { + id: "model.step.edit", + label: "Edit Step", + params: r#"{"name": "Bend", "procedure", "output", "boundary_conditions", "loads"} — as model.step.create; what is left out keeps its value"#, + run: step_edit, +}; + +pub const STEP_DELETE: CommandSpec = CommandSpec { + id: "model.step.delete", + label: "Delete Step", + params: r#"{"name": "Bend"}"#, + run: |s, p| s.model.remove_step(name(p)?).map(|_| Value::Null).map_err(Into::into), +}; + +pub const STEP_ACTIVATE: CommandSpec = CommandSpec { + id: "model.step.activate", + label: "Activate in Step", + params: r#"{"step": "Bend", "boundary_conditions": [names], "loads": [names], "active": true} — switches the listed items on (or off with false) in that step only"#, + run: step_activate, +}; + +pub const ACTIVITY: CommandSpec = CommandSpec { id: "model.activity", label: "Step Table", params: "{}", run: |s, _| Ok(activity(&s.model)) }; + +fn step_create(s: &mut Session, p: &Value) -> Result { + let n = name(p)?; + let mut step = Step::new(n, req(p, "procedure")?); + if let Some(o) = opt(p, "output")? { + step.output = o; + } + step.boundary_conditions = opt(p, "boundary_conditions")?.unwrap_or_default(); + step.loads = opt(p, "loads")?.unwrap_or_default(); + s.model.add_step(step, opt(p, "before")?)?; + Ok(json!({"name": n})) +} + +fn step_edit(s: &mut Session, p: &Value) -> Result { + let n = name(p)?; + let mut step = s.model.step(n).cloned().ok_or_else(|| EngineError::Params(format!("no step `{n}`")))?; + if let Some(procedure) = opt(p, "procedure")? { + step.procedure = procedure; + } + if let Some(o) = opt(p, "output")? { + step.output = o; + } + if let Some(b) = opt(p, "boundary_conditions")? { + step.boundary_conditions = b; + } + if let Some(l) = opt(p, "loads")? { + step.loads = l; + } + s.model.replace_step(step)?; + Ok(json!({"name": n})) +} + +fn step_activate(s: &mut Session, p: &Value) -> Result { + let step: &str = req(p, "step")?; + let active = opt(p, "active")?.unwrap_or(true); + let bcs: Vec<&str> = opt(p, "boundary_conditions")?.unwrap_or_default(); + let loads: Vec<&str> = opt(p, "loads")?.unwrap_or_default(); + if bcs.is_empty() && loads.is_empty() { + return Err(EngineError::Params("list `boundary_conditions` or `loads` to switch".into())); + } + for b in &bcs { + s.model.set_boundary_condition_active(step, b, active)?; + } + for l in &loads { + s.model.set_load_active(step, l, active)?; + } + Ok(activity(&s.model)) +} + +/// The step × boundary-condition/load table: the steps in order, and for each item the steps it +/// is active in. +fn activity(m: &Model) -> Value { + let row = + |active: &dyn Fn(&Step) -> bool| -> Value { m.steps().iter().filter(|s| active(s)).map(|s| s.name.as_str()).collect::>().into() }; + let bcs: serde_json::Map = + m.boundary_conditions().keys().map(|n| (n.clone(), row(&|s| s.boundary_conditions.contains(n)))).collect(); + let loads: serde_json::Map = m.loads().keys().map(|n| (n.clone(), row(&|s| s.loads.contains(n)))).collect(); + json!({ + "steps": m.steps().iter().map(|s| json!({"name": s.name, "procedure": s.procedure.name()})).collect::>(), + "boundary_conditions": bcs, + "loads": loads, + }) +} + +#[cfg(test)] +pub(crate) mod tests { + use serde_json::json; + + use crate::Session; + + /// A two-element bar tagged with body `Bar` and faces `left` and `right`. + pub(crate) fn tagged_bar(s: &mut Session) { + for (i, x) in [0.0, 1.0, 2.0].into_iter().enumerate() { + for (j, (y, z)) in [(0.0, 0.0), (0.0, 1.0), (1.0, 0.0), (1.0, 1.0)].into_iter().enumerate() { + s.execute("mesh.node.add", &json!({"label": 1 + i * 4 + j, "xyz": [x, y, z]})).unwrap(); + } + } + // Node (i, y, z) = 1 + 4i + 2y + z. + let n = |i: u32, y: u32, z: u32| 1 + 4 * i + 2 * y + z; + for e in 0..2 { + let nodes = [n(e, 0, 0), n(e + 1, 0, 0), n(e + 1, 1, 0), n(e, 1, 0), n(e, 0, 1), n(e + 1, 0, 1), n(e + 1, 1, 1), n(e, 1, 1)]; + s.execute("mesh.element.add", &json!({"label": e + 1, "kind": "hex8", "nodes": nodes})).unwrap(); + } + s.execute("mesh.tag", &json!({"geometry": {"body": "Bar"}, "elements": [1, 2]})).unwrap(); + s.execute("mesh.tag", &json!({"geometry": {"body": "Bar", "face": "left"}, "faces": [{"element": 1, "face": 6}]})).unwrap(); + s.execute("mesh.tag", &json!({"geometry": {"body": "Bar", "face": "right"}, "faces": [{"element": 2, "face": 4}]})).unwrap(); + } + + pub(crate) fn defined_bar(s: &mut Session) { + tagged_bar(s); + let script = json!([ + {"command": "model.units", "params": {"system": "mm-t-s-N-MPa"}}, + {"command": "model.set.create", "params": {"name": "all", "geometry": [{"body": "Bar"}]}}, + {"command": "model.set.create", "params": {"name": "root", "geometry": [{"body": "Bar", "face": "left"}]}}, + {"command": "model.set.create", "params": {"name": "end", "geometry": [{"body": "Bar", "face": "right"}]}}, + {"command": "model.material.create", "params": {"name": "Steel", "elastic": {"type": "isotropic", "young": 210000, "poisson": 0.3}, "density": 7.85e-9}}, + {"command": "model.section.create", "params": {"name": "bar", "set": "all", "material": "Steel", "kind": {"type": "solid"}, "formulation": {"incompatible_modes": true}}}, + {"command": "model.reference_point.create", "params": {"name": "RP", "xyz": [3, 0.5, 0.5]}}, + {"command": "model.coupling.create", "params": {"name": "pad", "reference_point": "RP", "set": "end", "kind": "kinematic"}}, + {"command": "model.boundary_condition.create", "params": {"name": "clamp", "type": "fixed", "target": {"set": "root"}}}, + {"command": "model.load.create", "params": {"name": "push", "type": "force", "target": {"reference_point": "RP"}, "vector": [0, 0, -100]}}, + {"command": "model.load.create", "params": {"name": "press", "type": "pressure", "set": "end", "magnitude": 2}}, + {"command": "model.step.create", "params": {"name": "Bend", "procedure": {"type": "static"}, "boundary_conditions": ["clamp"], "loads": ["push"]}}, + {"command": "model.step.create", "params": {"name": "Modes", "procedure": {"type": "frequency", "modes": 6}, "boundary_conditions": ["clamp"]}}, + ]); + s.run_script(&script).unwrap(); + } + + #[test] + fn a_script_defines_a_model_that_checks_clean() { + let mut s = Session::default(); + defined_bar(&mut s); + let check = s.execute("model.check", &json!({})).unwrap(); + assert_eq!(check["ready"], true, "{check}"); + // `press` is defined but in no step: a warning, not an error. + assert_eq!(check["warnings"], 1, "{check}"); + assert!(check["diagnostics"][0]["subject"] == "load `press`", "{check}"); + + let r = s.execute("model.set.resolve", &json!({"name": "end"})).unwrap(); + assert_eq!(r["nodes"], json!({"count": 4, "labels": [9, 10, 11, 12]})); + assert_eq!(r["faces"]["labels"], json!([{"element": 2, "face": 4}])); + + let table = s.execute("model.step.activate", &json!({"step": "Bend", "loads": ["press"]})).unwrap(); + assert_eq!(table["loads"]["press"], json!(["Bend"])); + assert_eq!(table["boundary_conditions"]["clamp"], json!(["Bend", "Modes"])); + assert_eq!(s.execute("model.check", &json!({})).unwrap()["warnings"], 0); + let off = s.execute("model.step.activate", &json!({"step": "Bend", "loads": ["press", "push"], "active": false})).unwrap(); + assert_eq!(off["loads"], json!({"press": [], "push": []})); + + let info = s.execute("model.info", &json!({})).unwrap(); + assert_eq!(info["units"], "mm-t-s-N-MPa"); + assert_eq!(info["definitions"]["sets"], 3); + assert_eq!(info["definitions"]["steps"], 2); + } + + #[test] + fn edits_and_deletes_keep_the_model_consistent() { + let mut s = Session::default(); + defined_bar(&mut s); + let journal = s.journal.len(); + + // Refused: a taken name, a missing reference, a referenced item, a typo in a field. + assert!(s.execute("model.set.create", &json!({"name": "all", "nodes": [1]})).is_err()); + assert!(s.execute("model.section.create", &json!({"name": "x", "set": "nope", "material": "Steel", "kind": {"type": "solid"}})).is_err()); + let e = s.execute("model.set.delete", &json!({"name": "end"})).unwrap_err().to_string(); + assert!(e.contains("coupling `pad`") && e.contains("load `press`"), "{e}"); + assert!( + s.execute("model.material.create", &json!({"name": "Al", "elastic": {"type": "isotropic", "young": 70000, "poison": 0.33}})).is_err() + ); + assert!(s.execute("model.step.edit", &json!({"name": "Bend", "loads": ["nope"]})).is_err()); + assert!(s.execute("model.step.create", &json!({"name": "Late", "procedure": {"type": "static"}, "before": "nope"})).is_err()); + assert_eq!(s.journal.len(), journal, "refusals change nothing"); + + // Edits replace; the step keeps its table entries unless told otherwise. + s.execute( + "model.material.edit", + &json!({"name": "Steel", "elastic": {"type": "isotropic", "young": 200000, "poisson": 0.29}, "density": 7.8e-9}), + ) + .unwrap(); + s.execute("model.step.edit", &json!({"name": "Modes", "procedure": {"type": "frequency", "modes": 12, "max_hz": 1000}})).unwrap(); + assert!(s.model.step("Modes").unwrap().boundary_conditions.contains("clamp")); + s.execute( + "model.step.create", + &json!({"name": "Pre", "procedure": {"type": "static"}, "before": "Bend", "boundary_conditions": ["clamp"], "loads": ["press"]}), + ) + .unwrap(); + assert_eq!(s.model.steps().iter().map(|st| st.name.as_str()).collect::>(), ["Pre", "Bend", "Modes"]); + + // A load listed by a step is not deleted until no step lists it. + assert!(s.execute("model.load.delete", &json!({"name": "press"})).unwrap_err().to_string().contains("step `Pre`")); + s.execute("model.step.activate", &json!({"step": "Pre", "loads": ["press"], "active": false})).unwrap(); + s.execute("model.load.delete", &json!({"name": "press"})).unwrap(); + s.execute("model.coupling.delete", &json!({"name": "pad"})).unwrap(); + s.execute("model.set.delete", &json!({"name": "end"})).unwrap(); + s.execute("model.step.delete", &json!({"name": "Pre"})).unwrap(); + assert!(s.execute("model.step.delete", &json!({"name": "Pre"})).is_err()); + } + + #[test] + fn the_check_names_the_bad_items_and_sets_follow_a_retag() { + let mut s = Session::default(); + defined_bar(&mut s); + // Point the section at a face set: it has no elements, and the elements lose their section. + s.execute("model.section.edit", &json!({"name": "bar", "set": "root", "material": "Steel", "kind": {"type": "solid"}})).unwrap(); + let check = s.execute("model.check", &json!({})).unwrap(); + assert_eq!(check["ready"], false); + let messages: Vec = check["diagnostics"] + .as_array() + .unwrap() + .iter() + .map(|d| format!("{}: {}", d["subject"].as_str().unwrap(), d["message"].as_str().unwrap())) + .collect(); + assert!(messages.contains(&"section `bar`: set `root` has no elements".to_string()), "{messages:#?}"); + assert!(messages.contains(&"model: elements 1, 2 have no section".to_string()), "{messages:#?}"); + + // Clearing the tags leaves geometry sets unresolved, and says so. + s.execute("mesh.tag.clear", &json!({})).unwrap(); + let e = s.execute("model.set.resolve", &json!({"name": "root"})).unwrap_err().to_string(); + assert!(e.contains("face left of body Bar is not on the mesh (the mesh carries no geometry tags)"), "{e}"); + } +} diff --git a/crates/engine/src/cmd/mesh.rs b/crates/engine/src/cmd/mesh.rs index 9265612..df4268a 100644 --- a/crates/engine/src/cmd/mesh.rs +++ b/crates/engine/src/cmd/mesh.rs @@ -7,7 +7,7 @@ use modelcraft_mesh::surface::{BoundarySurface, Defect, FeatureKind, SurfaceOpti use modelcraft_model::{ElementKind, MAX_ELEMENTS}; use serde_json::{Value, json}; -use super::{label, labels, optional_number, points, xyz}; +use super::{label, labels, listed, optional_number, points, xyz}; use crate::{CommandSpec, EngineError, Result, Session}; pub const NODE_ADD: CommandSpec = CommandSpec { @@ -96,18 +96,11 @@ fn element_remove(s: &mut Session, p: &Value) -> Result { Ok(Value::Null) } -/// Labels listed in a report: the first ones, and how many there are. -const LISTED: usize = 100; - /// A JSON number, or a string for ±∞ (JSON has no infinity). fn number(v: f64) -> Value { if v.is_finite() { json!(v) } else { json!(if v > 0.0 { "inf" } else { "-inf" }) } } -fn listed(v: &[u32]) -> Value { - json!({"count": v.len(), "labels": v.iter().take(LISTED).collect::>()}) -} - fn quality_report(s: &mut Session, p: &Value) -> Result { let bins = match p.get("bins") { None => 10, diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index 51ca70d..9a06fac 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -1,11 +1,14 @@ //! The command registry. Each area has its own module; [`ALL`] lists every command. +pub(crate) mod define; mod geometry; mod mesh; mod model; +mod tags; use modelcraft_model::MAX_LABEL; -use serde_json::Value; +use serde::Deserialize; +use serde_json::{Value, json}; use crate::{CommandSpec, EngineError, Result}; @@ -22,6 +25,37 @@ pub static ALL: &[CommandSpec] = &[ mesh::QUALITY, mesh::DELAUNAY_POINTS, mesh::SURFACE_CHECK, + tags::TAG, + tags::TAG_CLEAR, + model::UNITS, + model::CHECK, + model::SET_RESOLVE, + define::SET_CREATE, + define::SET_EDIT, + define::SET_DELETE, + define::MATERIAL_CREATE, + define::MATERIAL_EDIT, + define::MATERIAL_DELETE, + define::SECTION_CREATE, + define::SECTION_EDIT, + define::SECTION_DELETE, + define::REFERENCE_POINT_CREATE, + define::REFERENCE_POINT_EDIT, + define::REFERENCE_POINT_DELETE, + define::COUPLING_CREATE, + define::COUPLING_EDIT, + define::COUPLING_DELETE, + define::BOUNDARY_CONDITION_CREATE, + define::BOUNDARY_CONDITION_EDIT, + define::BOUNDARY_CONDITION_DELETE, + define::LOAD_CREATE, + define::LOAD_EDIT, + define::LOAD_DELETE, + define::STEP_CREATE, + define::STEP_EDIT, + define::STEP_DELETE, + define::STEP_ACTIVATE, + define::ACTIVITY, ]; /// A required label parameter (1..=u32::MAX; the model checks its own range). @@ -72,3 +106,25 @@ pub(crate) fn optional_number(p: &Value, key: &str, default: f64) -> Result Some(v) => v.as_f64().ok_or_else(|| EngineError::Params(format!("`{key}` must be a number"))), } } + +/// The `name` parameter. +pub(crate) fn name(p: &Value) -> Result<&str> { + p.get("name").and_then(Value::as_str).ok_or_else(|| EngineError::Params("`name` must be a name".into())) +} + +/// An optional parameter read as `T`'s JSON form; `None` when it is absent. +pub(crate) fn opt<'a, T: Deserialize<'a>>(p: &'a Value, key: &str) -> Result> { + p.get(key).map(|v| T::deserialize(v).map_err(|e| EngineError::Params(format!("`{key}`: {e}")))).transpose() +} + +/// A required parameter read as `T`'s JSON form. +pub(crate) fn req<'a, T: Deserialize<'a>>(p: &'a Value, key: &str) -> Result { + opt(p, key)?.ok_or_else(|| EngineError::Params(format!("`{key}` is required"))) +} + +/// Entries listed in a report: the first ones, and how many there are. +pub(crate) const LISTED: usize = 100; + +pub(crate) fn listed(v: &[T]) -> Value { + json!({"count": v.len(), "labels": v.get(..LISTED).unwrap_or(v)}) +} diff --git a/crates/engine/src/cmd/model.rs b/crates/engine/src/cmd/model.rs index 347ef9f..422a8cb 100644 --- a/crates/engine/src/cmd/model.rs +++ b/crates/engine/src/cmd/model.rs @@ -1,10 +1,62 @@ +use modelcraft_model::{Severity, UnitSystem}; use serde_json::{Value, json}; -use crate::{CommandSpec, Result, Session}; +use super::{listed, name}; +use crate::{CommandSpec, EngineError, Result, Session}; pub const INFO: CommandSpec = CommandSpec { id: "model.info", label: "Model Info", params: "{}", run: info }; +pub const UNITS: CommandSpec = CommandSpec { + id: "model.units", + label: "Unit System", + params: r#"{"system": "m-kg-s-N-Pa|mm-t-s-N-MPa|mm-kg-ms-kN-GPa|mm-g-ms-N-MPa|in-lbf-s-psi|unspecified"} — labels the system the model's numbers are in; nothing is converted"#, + run: units, +}; + +pub const CHECK: CommandSpec = CommandSpec { id: "model.check", label: "Check Model", params: "{}", run: check }; + +pub const SET_RESOLVE: CommandSpec = CommandSpec { + id: "model.set.resolve", + label: "Resolve Set", + params: r#"{"name": "root"} — the set's nodes, elements and element faces on the current mesh (the first 100 of each); fails naming every member that is not on the mesh"#, + run: set_resolve, +}; + fn info(s: &mut Session, _: &Value) -> Result { let counts: serde_json::Map = s.model.element_counts().into_iter().map(|(k, n)| (k.name().to_string(), json!(n))).collect(); - Ok(json!({"nodes": s.model.nodes().len(), "elements": counts})) + let definitions: serde_json::Map = s.model.definition_counts().into_iter().map(|(k, n)| (k.to_string(), json!(n))).collect(); + Ok(json!({"nodes": s.model.nodes().len(), "elements": counts, "units": s.model.units().label(), "definitions": definitions})) +} + +fn units(s: &mut Session, p: &Value) -> Result { + let label = p.get("system").and_then(Value::as_str).unwrap_or_default(); + let system = UnitSystem::from_label(label).ok_or_else(|| { + let known: Vec<&str> = UnitSystem::ALL.iter().map(|u| u.label()).collect(); + EngineError::Params(format!("`system` must be one of {}", known.join(", "))) + })?; + s.model.set_units(system); + // SolidCraft designs are in millimetres. + let warning = match system.length() { + Some(length) if length != "mm" && s.design.is_some() => { + Some(format!("the open design is in millimetres but this system's length unit is {length}")) + } + _ => None, + }; + Ok(json!({"system": system.label(), "warning": warning})) +} + +fn check(s: &mut Session, _: &Value) -> Result { + let diagnostics = s.model.check(); + let errors = diagnostics.iter().filter(|d| d.severity == Severity::Error).count(); + Ok(json!({ + "ready": errors == 0, + "errors": errors, + "warnings": diagnostics.len() - errors, + "diagnostics": diagnostics, + })) +} + +fn set_resolve(s: &mut Session, p: &Value) -> Result { + let r = s.model.resolve_set(name(p)?)?; + Ok(json!({"nodes": listed(&r.nodes), "elements": listed(&r.elements), "faces": listed(&r.faces)})) } diff --git a/crates/engine/src/cmd/tags.rs b/crates/engine/src/cmd/tags.rs new file mode 100644 index 0000000..156f2f9 --- /dev/null +++ b/crates/engine/src/cmd/tags.rs @@ -0,0 +1,63 @@ +//! Geometry tags on the mesh: which body, face, edge or vertex each mesh entity lies on. Meshers +//! write them; these commands set them by hand (for a mesh from elsewhere) and clear them. + +use modelcraft_model::{ElementFace, GeometryRef, MAX_ELEMENTS}; +use serde_json::{Value, json}; + +use super::{opt, req}; +use crate::{CommandSpec, EngineError, Result, Session}; + +pub const TAG: CommandSpec = CommandSpec { + id: "mesh.tag", + label: "Tag Mesh with Geometry", + params: r#"{"geometry": {"body": "Body1"} | {"body", "face"} | {"body", "edge"} | {"body", "vertex"}, "elements": [labels] (filling a body, face or edge), "faces": [{"element": 1, "face": 3}] (on a face), "edges": [{"element": 1, "edge": 2}] (on an edge), "nodes": [labels] (on a vertex)} — records where the mesh lies on the geometry, so sets named by geometry resolve"#, + run: tag, +}; + +pub const TAG_CLEAR: CommandSpec = CommandSpec { id: "mesh.tag.clear", label: "Clear Geometry Tags", params: "{}", run: tag_clear }; + +#[derive(serde::Deserialize)] +#[serde(deny_unknown_fields)] +struct ElementEdge { + element: u32, + edge: u8, +} + +/// An optional list parameter of at most [`MAX_ELEMENTS`] entries. +fn list<'a, T: serde::Deserialize<'a>>(p: &'a Value, key: &str) -> Result> { + if p.get(key).and_then(Value::as_array).is_some_and(|a| a.len() > MAX_ELEMENTS) { + return Err(EngineError::Params(format!("`{key}` must be a list of at most {MAX_ELEMENTS} entries"))); + } + Ok(opt(p, key)?.unwrap_or_default()) +} + +fn tag(s: &mut Session, p: &Value) -> Result { + let geometry: GeometryRef = req(p, "geometry")?; + let elements: Vec = list(p, "elements")?; + let nodes: Vec = list(p, "nodes")?; + let faces: Vec = list(p, "faces")?; + let edges: Vec = list(p, "edges")?; + let count = elements.len() + nodes.len() + faces.len() + edges.len(); + if count == 0 { + return Err(EngineError::Params("list `elements`, `faces`, `edges` or `nodes` to tag".into())); + } + let t = s.model.geometry_tag(&geometry)?; + for e in elements { + s.model.tag_element(e, t)?; + } + for f in faces { + s.model.tag_element_face(f.element, f.face, t)?; + } + for e in edges { + s.model.tag_element_edge(e.element, e.edge, t)?; + } + for n in nodes { + s.model.tag_node(n, t)?; + } + Ok(json!({"geometry": geometry, "tagged": count})) +} + +fn tag_clear(s: &mut Session, _: &Value) -> Result { + s.model.clear_geometry_tags(); + Ok(Value::Null) +} diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 8b5d2e6..07babef 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -188,6 +188,41 @@ mod tests { json!({"surface": -1}), json!({"surface": 1e300}), json!({"surface": "0", "weld": -1, "feature_angle": 1e300}), + // Model definitions. + json!({"name": "x".repeat(100_000)}), + json!({"name": "A\u{0}"}), + json!({"name": 7, "nodes": [1]}), + json!({"name": "A", "nodes": [0, -1, 1e300]}), + json!({"name": "A", "nodes": vec![u32::MAX; 3], "faces": [{"element": 1, "face": 0}]}), + json!({"name": "A", "geometry": [{"body": ""}, {"body": "B", "face": "f", "edge": "e"}]}), + json!({"name": "A", "geometry": [{"body": "x".repeat(10_000)}]}), + json!({"name": "A", "elastic": {"type": "isotropic", "young": 1e308, "poisson": 0.49999999}}), + json!({"name": "A", "elastic": {"type": "engineering_constants", "young": [1e-300, 1e300, 1], "poisson": [1e300, 0, 0], "shear": [0, 0, 0]}}), + json!({"name": "A", "elastic": {"type": "isotropic", "young": "stiff"}, "density": -0.0}), + json!({"name": "A", "set": "A", "material": "A", "kind": {"type": "solid"}}), + json!({"name": "A", "set": "all", "material": "Steel", "kind": {"type": "beam", "profile": {"type": "pipe", "radius": 1, "thickness": 2}, "orientation": [0, 0, 0]}}), + json!({"name": "A", "set": "all", "material": "Steel", "kind": {"type": "shell", "thickness": 1e308, "offset": -0.5, "points": 4294967295u64}}), + json!({"name": "A", "xyz": [1e308, 0, 0]}), + json!({"name": "A", "reference_point": "RP", "set": "end", "kind": "rigid", "dofs": [0, 7]}), + json!({"name": "A", "type": "fixed", "target": {"set": "nope"}}), + json!({"name": "A", "type": "displacement", "target": {"reference_point": "RP"}, "values": [null, null, null, null, null, null]}), + json!({"name": "A", "type": "displacement", "target": {"set": "root"}, "values": [1e308, -1e308, 0, 0, 0, 0, 0]}), + json!({"name": "A", "type": "pressure", "set": "end", "magnitude": -1e308}), + json!({"name": "A", "type": "gravity", "acceleration": [0, 0, 0]}), + json!({"name": "A", "type": "force", "target": {"set": "end"}, "vector": [1, 0, 0], "distribution": "lumpy"}), + json!({"name": "A", "type": "traction", "set": "root", "vector": [1e308, 1e308, 1e308]}), + json!({"name": "A", "type": "traction", "set": "nope", "vector": [1, 0, 0]}), + json!({"name": "RP", "xyz": [0, 0, 0], "label": 4294967295u64}), + json!({"name": "A", "xyz": [0, 0, 0], "label": 1}), + json!({"name": "A", "type": "body_force", "set": "all", "vector": [1, 2]}), + json!({"name": "A", "procedure": {"type": "frequency", "modes": 4294967295u64, "min_hz": -1}}), + json!({"name": "A", "procedure": {"type": "frequency", "modes": 5}, "output": ["everything"], "before": 3}), + json!({"name": "Bend", "procedure": {"type": "static"}, "loads": "push"}), + json!({"step": "Bend", "loads": ["nope"], "active": "yes"}), + json!({"step": "Bend", "boundary_conditions": [], "loads": []}), + json!({"system": "furlongs"}), + json!({"geometry": {"body": "Bar"}, "elements": [1, 99], "faces": [{"element": 1, "face": 200}]}), + json!({"geometry": {"body": "Bar", "vertex": "v"}, "nodes": [1], "edges": [{"element": 1, "edge": 0}]}), ]; for spec in commands() { for p in &hostile { @@ -197,6 +232,19 @@ mod tests { } assert!(Session::default().execute("no.such", &json!({})).is_err()); + // The same against a defined model, so edits and deletes reach existing definitions. + let mut defined = Session::default(); + crate::cmd::define::tests::defined_bar(&mut defined); + for spec in commands() { + for p in &hostile { + let mut s = defined.clone(); + let before = s.model.clone(); + if s.execute(spec.id, p).is_err() { + assert_eq!(s.model, before, "{} {p}", spec.id); + } + } + } + // The same against an open design. let (_dir, path) = sample_design_file(&Plate::default()); let mut opened = Session::default(); diff --git a/crates/model/src/bcs.rs b/crates/model/src/bcs.rs new file mode 100644 index 0000000..9d13694 --- /dev/null +++ b/crates/model/src/bcs.rs @@ -0,0 +1,238 @@ +//! Boundary conditions: prescribed degrees of freedom on a set's nodes or a reference point. +//! +//! Degrees of freedom are numbered 1–3 for translations along X, Y, Z and 4–6 for rotations about +//! them, in the global axes. Which steps a boundary condition is active in is listed explicitly by +//! each step ([`crate::steps`]). + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::{Error, Model, Result, context, finite, insert_new, not_in_use, remove_existing, replace_existing, require}; + +/// A set of degrees of freedom (1–6), written as a sorted list such as `[1, 2, 3]`. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(try_from = "Vec", into = "Vec")] +pub struct Dofs(u8); + +impl Dofs { + pub fn all() -> Dofs { + Dofs(0b111_111) + } + + pub fn none() -> Dofs { + Dofs(0) + } + + pub fn contains(self, dof: u8) -> bool { + (1..=6).contains(&dof) && self.0 & (1 << (dof - 1)) != 0 + } + + pub fn is_empty(self) -> bool { + self.0 == 0 + } + + /// The degrees of freedom in order. + pub fn list(self) -> Vec { + (1..=6).filter(|d| self.contains(*d)).collect() + } +} + +impl TryFrom> for Dofs { + type Error = Error; + fn try_from(v: Vec) -> Result { + let mut mask = 0u8; + for d in v { + if !(1..=6).contains(&d) { + return Err(Error::Invalid(format!("degree of freedom {d} is not in 1..=6"))); + } + mask |= 1 << (d - 1); + } + Ok(Dofs(mask)) + } +} + +impl From for Vec { + fn from(d: Dofs) -> Vec { + d.list() + } +} + +/// What a boundary condition, force or moment acts on. +#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum Target { + /// The nodes of a set. + Set(String), + /// A reference point. + ReferencePoint(String), +} + +impl Target { + pub fn set(&self) -> Option<&str> { + match self { + Target::Set(s) => Some(s), + Target::ReferencePoint(_) => None, + } + } + + pub fn reference_point(&self) -> Option<&str> { + match self { + Target::ReferencePoint(r) => Some(r), + Target::Set(_) => None, + } + } + + pub(crate) fn check_exists(&self, model: &Model) -> Result<()> { + match self { + Target::Set(s) => require(&model.sets, "set", s), + Target::ReferencePoint(r) => require(&model.reference_points, "reference point", r), + } + } +} + +impl std::fmt::Display for Target { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Target::Set(s) => write!(f, "set `{s}`"), + Target::ReferencePoint(r) => write!(f, "reference point `{r}`"), + } + } +} + +/// The plane of a symmetry boundary condition, named by its normal axis. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum Symmetry { + /// The plane normal to X: no X translation, no rotation about Y or Z. + X, + Y, + Z, +} + +impl Symmetry { + /// The fixed degrees of freedom. + pub fn dofs(self) -> Dofs { + match self { + Symmetry::X => Dofs(0b110_001), + Symmetry::Y => Dofs(0b101_010), + Symmetry::Z => Dofs(0b011_100), + } + } +} + +/// A boundary condition. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] +pub enum BoundaryCondition { + /// Every degree of freedom held at zero. + Fixed { target: Target }, + /// Prescribed values for some degrees of freedom: `values[i]` is degree of freedom `i + 1` + /// (`null` leaves it free, `0` holds it). + Displacement { target: Target, values: [Option; 6] }, + /// Mirror symmetry about a plane through the nodes. + Symmetry { target: Target, plane: Symmetry }, +} + +impl BoundaryCondition { + pub fn target(&self) -> &Target { + match self { + BoundaryCondition::Fixed { target } | BoundaryCondition::Displacement { target, .. } | BoundaryCondition::Symmetry { target, .. } => { + target + } + } + } + + /// The constrained degrees of freedom. + pub fn dofs(&self) -> Dofs { + match self { + BoundaryCondition::Fixed { .. } => Dofs::all(), + BoundaryCondition::Displacement { values, .. } => { + Dofs(values.iter().enumerate().filter(|(_, v)| v.is_some()).fold(0, |m, (i, _)| m | (1 << i))) + } + BoundaryCondition::Symmetry { plane, .. } => plane.dofs(), + } + } + + /// The prescribed value of a constrained degree of freedom (1–6); `None` when it is free. + pub fn value(&self, dof: u8) -> Option { + match self { + BoundaryCondition::Displacement { values, .. } => values.get(usize::from(dof).checked_sub(1)?).copied().flatten(), + _ => self.dofs().contains(dof).then_some(0.0), + } + } + + pub fn validate(&self) -> Result<()> { + if let BoundaryCondition::Displacement { values, .. } = self { + for (i, v) in values.iter().enumerate() { + if let Some(v) = v { + finite(&format!("degree of freedom {}", i + 1), *v)?; + } + } + if values.iter().all(Option::is_none) { + return Err(Error::Invalid("a displacement boundary condition needs at least one value".into())); + } + } + Ok(()) + } +} + +impl Model { + pub fn boundary_conditions(&self) -> &BTreeMap { + &self.boundary_conditions + } + + pub fn boundary_condition(&self, name: &str) -> Option<&BoundaryCondition> { + self.boundary_conditions.get(name) + } + + fn check_bc(&self, name: &str, bc: &BoundaryCondition) -> Result<()> { + bc.validate().map_err(|e| context("boundary condition", name, e))?; + bc.target().check_exists(self).map_err(|e| context("boundary condition", name, e)) + } + + pub fn add_boundary_condition(&mut self, name: &str, bc: BoundaryCondition) -> Result<()> { + self.check_bc(name, &bc)?; + insert_new(&mut self.boundary_conditions, "boundary condition", name, bc) + } + + pub fn replace_boundary_condition(&mut self, name: &str, bc: BoundaryCondition) -> Result<()> { + self.check_bc(name, &bc)?; + replace_existing(&mut self.boundary_conditions, "boundary condition", name, bc) + } + + /// Remove a boundary condition that no step lists. + pub fn remove_boundary_condition(&mut self, name: &str) -> Result<()> { + let users = self.steps.iter().filter(|s| s.boundary_conditions.contains(name)).map(|s| format!("step `{}`", s.name)).collect(); + not_in_use("boundary condition", name, users)?; + remove_existing(&mut self.boundary_conditions, "boundary condition", name).map(|_| ()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn dofs_and_values() { + let d: Dofs = serde_json::from_value(serde_json::json!([3, 1, 1])).unwrap(); + assert_eq!(d.list(), vec![1, 3]); + assert_eq!(serde_json::to_value(d).unwrap(), serde_json::json!([1, 3])); + assert!(serde_json::from_value::(serde_json::json!([0])).is_err()); + assert!(serde_json::from_value::(serde_json::json!([7])).is_err()); + assert_eq!(Symmetry::X.dofs().list(), vec![1, 5, 6]); + assert_eq!(Symmetry::Y.dofs().list(), vec![2, 4, 6]); + assert_eq!(Symmetry::Z.dofs().list(), vec![3, 4, 5]); + + let t = Target::Set("s".into()); + let bc = BoundaryCondition::Displacement { target: t.clone(), values: [Some(0.0), None, Some(1.5), None, None, None] }; + assert_eq!(bc.dofs().list(), vec![1, 3]); + assert_eq!((bc.value(1), bc.value(2), bc.value(3), bc.value(0), bc.value(9)), (Some(0.0), None, Some(1.5), None, None)); + assert_eq!(BoundaryCondition::Fixed { target: t.clone() }.value(6), Some(0.0)); + assert!(BoundaryCondition::Displacement { target: t.clone(), values: [None; 6] }.validate().is_err()); + assert!(BoundaryCondition::Displacement { target: t, values: [Some(f64::NAN), None, None, None, None, None] }.validate().is_err()); + let json = serde_json::json!({"type": "symmetry", "target": {"reference_point": "RP"}, "plane": "z"}); + let parsed: BoundaryCondition = serde_json::from_value(json.clone()).unwrap(); + assert_eq!(serde_json::to_value(&parsed).unwrap(), json); + } +} diff --git a/crates/model/src/check.rs b/crates/model/src/check.rs new file mode 100644 index 0000000..aedd06a --- /dev/null +++ b/crates/model/src/check.rs @@ -0,0 +1,436 @@ +//! The model check and the resolved view. +//! +//! [`Model::check`] lists everything that keeps the model from being analysed (errors) or looks +//! wrong (warnings), each naming the item it is about. [`Model::resolve`] gives the view the deck +//! writers and the solver read: every set resolved against the current mesh, the section of every +//! element and the node label of every reference point; it fails when the check finds an error. + +use std::collections::{BTreeMap, BTreeSet}; + +use serde::{Deserialize, Serialize}; + +use crate::loads::Distribution; +use crate::{BoundaryCondition, Elastic, Error, Load, Model, Procedure, Result, SectionKind, Target}; + +/// How serious a [`Diagnostic`] is. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum Severity { + /// The model cannot be analysed until it is fixed. + Error, + /// The model can be analysed but something looks wrong. + Warning, +} + +/// One finding of the check. +#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Serialize, Deserialize)] +pub struct Diagnostic { + pub severity: Severity, + /// The item it is about, such as ``section `web` `` or `model`. + pub subject: String, + pub message: String, +} + +impl std::fmt::Display for Diagnostic { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let s = match self.severity { + Severity::Error => "error", + Severity::Warning => "warning", + }; + write!(f, "{s}: {}: {}", self.subject, self.message) + } +} + +/// A set resolved against the mesh. Every list is sorted and has no repeats. +#[derive(Debug, Clone, Default, PartialEq, Eq, Serialize, Deserialize)] +pub struct ResolvedSet { + pub nodes: Vec, + pub elements: Vec, + pub faces: Vec, +} + +impl ResolvedSet { + pub fn is_empty(&self) -> bool { + self.nodes.is_empty() && self.elements.is_empty() && self.faces.is_empty() + } +} + +/// The model resolved against its mesh: what writers and the solver read besides the +/// definitions themselves. When resolution succeeds, every element is in exactly one section's +/// set (`sets[§ion.set].elements`). +#[derive(Debug, Clone, Default, PartialEq)] +pub struct Resolved { + /// Every set, by name. + pub sets: BTreeMap, + /// The node label of every reference point ([`Model::reference_point_nodes`]). + pub reference_point_nodes: BTreeMap, + /// Forces, moments and tractions as nodal values, sorted by node: what writers and the solver + /// apply, so every one of them distributes a load the same way. A force or traction gives + /// force components, a moment moment components. + pub nodal_loads: BTreeMap>, + /// The check's warnings (a resolved model has no errors). + pub warnings: Vec, +} + +impl Resolved { + /// The nodes a boundary condition, force or moment acts on. + pub fn target_nodes(&self, target: &Target) -> &[u32] { + match target { + Target::Set(s) => self.sets.get(s).map(|r| r.nodes.as_slice()).unwrap_or_default(), + Target::ReferencePoint(r) => self.reference_point_nodes.get(r).map(std::slice::from_ref).unwrap_or_default(), + } + } +} + +/// Labels listed in a message before it says how many more there are. +const LISTED: usize = 5; + +fn first_labels(labels: &[u32]) -> String { + let shown: Vec = labels.iter().take(LISTED).map(u32::to_string).collect(); + match labels.len().saturating_sub(LISTED) { + 0 => shown.join(", "), + more => format!("{}, and {more} more", shown.join(", ")), + } +} + +#[derive(Default)] +struct Diagnostics(Vec); + +impl Diagnostics { + fn error(&mut self, subject: impl Into, message: impl Into) { + self.0.push(Diagnostic { severity: Severity::Error, subject: subject.into(), message: message.into() }); + } + + fn warning(&mut self, subject: impl Into, message: impl Into) { + self.0.push(Diagnostic { severity: Severity::Warning, subject: subject.into(), message: message.into() }); + } + + /// A warning when a magnitude check found one. + fn warning_if(&mut self, subject: &str, message: Option) { + if let Some(m) = message { + self.warning(subject, m); + } + } +} + +impl Model { + /// Everything that keeps the model from being analysed or looks wrong, errors first, each + /// naming its item. + pub fn check(&self) -> Vec { + self.analyze().1 + } + + /// Resolve the model for writing or solving. Fails with every error the check finds. + pub fn resolve(&self) -> Result { + let (mut resolved, diagnostics) = self.analyze(); + let errors: Vec = diagnostics.iter().filter(|d| d.severity == Severity::Error).map(ToString::to_string).collect(); + if !errors.is_empty() { + return Err(Error::Invalid(format!("the model is not ready for analysis: {}", errors.join("; ")))); + } + resolved.warnings = diagnostics; + Ok(resolved) + } + + fn analyze(&self) -> (Resolved, Vec) { + let mut d = Diagnostics::default(); + let mut resolved = Resolved::default(); + let units = self.units; + + if self.elements.is_empty() { + d.error("model", "the mesh has no elements"); + } + if units.length().is_none() { + d.warning("model", "no unit system is chosen, so magnitudes are not checked"); + } + + // Sets. + for (name, set) in &self.sets { + let (r, missing) = self.resolve_definition(set); + for m in missing { + d.error(format!("set `{name}`"), m); + } + if r.is_empty() { + if self.set_users(name).is_empty() { + d.warning(format!("set `{name}`"), "resolves to nothing on the current mesh"); + } else { + d.error(format!("set `{name}`"), "resolves to nothing on the current mesh"); + } + } + resolved.sets.insert(name.clone(), r); + } + let empty = crate::ResolvedSet::default(); + let set_of = |name: &str| resolved.sets.get(name).unwrap_or(&empty); + + // Materials: magnitudes for the unit system. + for (name, m) in &self.materials { + let subject = format!("material `{name}`"); + match m.elastic { + Some(Elastic::Isotropic { young, .. }) => d.warning_if(&subject, units.check_modulus("Young's modulus", young)), + Some(Elastic::EngineeringConstants { young, shear, .. }) => { + for (i, e) in young.iter().enumerate() { + d.warning_if(&subject, units.check_modulus(&format!("E{}", i + 1), *e)); + } + for (n, g) in crate::materials::SHEAR_MODULI.iter().zip(shear) { + d.warning_if(&subject, units.check_modulus(n, g)); + } + } + None => {} + } + if let Some(rho) = m.density { + d.warning_if(&subject, units.check_density(rho)); + } + } + + // Sections: element kinds, formulation, one section per element. + let mut owner: BTreeMap = BTreeMap::new(); + for (name, sec) in &self.sections { + let subject = format!("section `{name}`"); + let r = set_of(&sec.set); + if r.elements.is_empty() { + d.error(&subject, format!("set `{}` has no elements", sec.set)); + } + let fits = |k: crate::ElementKind| match sec.kind { + SectionKind::Solid {} => k.dimension() == 3, + SectionKind::Shell { .. } => k.dimension() == 2, + SectionKind::Beam { .. } => k == crate::ElementKind::Beam2, + SectionKind::Rod { .. } => k == crate::ElementKind::Rod2, + }; + let mut wrong = Vec::new(); + let mut not_hex8 = Vec::new(); + let mut along = Vec::new(); + let mut taken: BTreeMap<&str, Vec> = BTreeMap::new(); + for &e in &r.elements { + let Some(el) = self.elements.get(&e) else { continue }; + if !fits(el.kind) { + wrong.push(e); + } + if sec.formulation.incompatible_modes && el.kind != crate::ElementKind::Hex8 { + not_hex8.push(e); + } + if let SectionKind::Beam { orientation, .. } = sec.kind + && self.parallel_to_axis(&el.nodes, orientation) + { + along.push(e); + } + match owner.get(&e) { + Some(other) => taken.entry(other).or_default().push(e), + None => { + owner.insert(e, name); + } + } + } + if !wrong.is_empty() { + d.error(&subject, format!("a {} section cannot hold elements {}", sec.kind.name(), first_labels(&wrong))); + } + if !not_hex8.is_empty() { + d.error(&subject, format!("incompatible modes apply to hex8 elements only, not elements {}", first_labels(¬_hex8))); + } + if !along.is_empty() { + d.error(&subject, format!("the orientation vector is parallel to elements {}", first_labels(&along))); + } + for (other, labels) in taken { + d.error(&subject, format!("elements {} are also in section `{other}`", first_labels(&labels))); + } + if self.materials.get(&sec.material).is_some_and(|m| m.elastic.is_none()) { + d.error(&subject, format!("material `{}` has no elastic constants", sec.material)); + } + } + let unsectioned: Vec = self.elements.keys().filter(|e| !owner.contains_key(e)).copied().collect(); + if !unsectioned.is_empty() { + d.error("model", format!("elements {} have no section", first_labels(&unsectioned))); + } + + // Mass: frequency steps and gravity need a density in every material in use. + let needs_mass = self + .steps + .iter() + .find(|s| matches!(s.procedure, Procedure::Frequency { .. })) + .map(|s| format!("step `{}` extracts frequencies", s.name)) + .or_else(|| self.loads.iter().find(|(_, l)| matches!(l, Load::Gravity { .. })).map(|(n, _)| format!("load `{n}` is gravity"))); + if let Some(why) = needs_mass { + let used: BTreeSet<&str> = self.sections.values().map(|s| s.material.as_str()).collect(); + for m in used { + if self.materials.get(m).is_some_and(|m| m.density.is_none()) { + d.error(format!("material `{m}`"), format!("has no density, and {why}")); + } + } + } + + // Reference points and couplings. + resolved.reference_point_nodes = self.reference_point_nodes(); + for (name, label) in &resolved.reference_point_nodes { + if self.nodes.contains_key(label) { + d.error(format!("reference point `{name}`"), format!("its node label {label} is also a mesh node")); + } + } + for (name, c) in &self.couplings { + if set_of(&c.set).nodes.is_empty() { + d.error(format!("coupling `{name}`"), format!("set `{}` has no nodes", c.set)); + } + } + + // Nodes with rotational degrees of freedom: those of shells and beams, and reference points + // (worked out only when a boundary condition or load needs it). + let rotational = std::cell::OnceCell::new(); + let has_rotations = |t: &Target| match t { + Target::ReferencePoint(_) => true, + Target::Set(s) => { + let rotational: &BTreeSet = rotational.get_or_init(|| { + self.elements + .values() + .filter(|e| e.kind.dimension() == 2 || e.kind == crate::ElementKind::Beam2) + .flat_map(|e| e.nodes.iter().copied()) + .collect() + }); + set_of(s).nodes.iter().any(|n| rotational.contains(n)) + } + }; + + // Boundary conditions. + for (name, bc) in &self.boundary_conditions { + let subject = format!("boundary condition `{name}`"); + if let Target::Set(s) = bc.target() + && set_of(s).nodes.is_empty() + { + d.error(&subject, format!("set `{s}` has no nodes")); + } + if let BoundaryCondition::Displacement { values, target } = bc + && values.get(3..).is_some_and(|r| r.iter().any(Option::is_some)) + && !has_rotations(target) + { + d.warning(&subject, format!("{target} has only solid-element nodes, which have no rotations to prescribe")); + } + } + + // Loads. + let mut nodal_loads = BTreeMap::new(); + for (name, load) in &self.loads { + let subject = format!("load `{name}`"); + match self.nodal_values(load, &resolved) { + Ok(Some(values)) => { + nodal_loads.insert(name.clone(), values); + } + Ok(None) => {} + Err(m) => d.error(&subject, m), + } + match load { + Load::Moment { target, .. } if !has_rotations(target) => { + d.error(&subject, format!("{target} has no nodes with rotations; apply the moment to a reference point coupled to it")) + } + Load::Pressure { set, .. } if set_of(set).faces.is_empty() => { + d.error(&subject, format!("set `{set}` has no element faces to press on")) + } + Load::BodyForce { set, .. } if set_of(set).elements.is_empty() => d.error(&subject, format!("set `{set}` has no elements")), + Load::Gravity { acceleration, set } => { + let g = acceleration.iter().map(|c| c * c).sum::().sqrt(); + d.warning_if(&subject, units.check_gravity(g)); + if let Some(s) = set + && set_of(s).elements.is_empty() + { + d.error(&subject, format!("set `{s}` has no elements")); + } + } + _ => {} + } + } + + // Steps and the activity table. + if self.steps.is_empty() { + d.error("model", "there is no analysis step"); + } + for s in &self.steps { + let subject = format!("step `{}`", s.name); + match s.procedure { + Procedure::Static {} => { + let prescribes_motion = s + .boundary_conditions + .iter() + .filter_map(|b| self.boundary_conditions.get(b)) + .any(|b| (1..=6).any(|dof| b.value(dof).is_some_and(|v| v != 0.0))); + if s.loads.is_empty() && !prescribes_motion { + d.warning(&subject, "nothing loads the model in this step"); + } + } + Procedure::Frequency { .. } => { + if !s.loads.is_empty() { + d.warning(&subject, "loads have no effect on natural frequencies"); + } + } + } + if s.boundary_conditions.is_empty() { + d.warning(&subject, "no boundary condition is active, so the model is free to move as a rigid body"); + } + } + let active_bcs: BTreeSet<&String> = self.steps.iter().flat_map(|s| &s.boundary_conditions).collect(); + let active_loads: BTreeSet<&String> = self.steps.iter().flat_map(|s| &s.loads).collect(); + for n in self.boundary_conditions.keys().filter(|n| !active_bcs.contains(n)) { + d.warning(format!("boundary condition `{n}`"), "is not active in any step"); + } + for n in self.loads.keys().filter(|n| !active_loads.contains(n)) { + d.warning(format!("load `{n}`"), "is not active in any step"); + } + + resolved.nodal_loads = nodal_loads; + let mut out = d.0; + out.sort_by_key(|a| a.severity); + (resolved, out) + } + + /// A force, moment or traction as nodal values (`None` for other loads), or why it cannot be + /// distributed. + fn nodal_values(&self, load: &Load, r: &Resolved) -> std::result::Result>, String> { + let scale = |v: [f64; 3], k: f64| v.map(|c| c * k); + let (target, vector, distribution) = match load { + Load::Force { target, vector, distribution } | Load::Moment { target, vector, distribution } => (target, *vector, *distribution), + Load::Traction { set, vector } => { + let faces = r.sets.get(set).map(|s| s.faces.as_slice()).unwrap_or_default(); + if faces.is_empty() { + return Err(format!("set `{set}` has no element faces to load")); + } + let w = self.surface_weights(faces).map_err(|e| e.to_string())?; + return Ok(Some(w.into_iter().map(|(n, a)| (n, scale(*vector, a))).collect())); + } + _ => return Ok(None), + }; + let empty = ResolvedSet::default(); + let set = match target { + Target::ReferencePoint(rp) => { + return r + .reference_point_nodes + .get(rp) + .map(|n| Some(vec![(*n, vector)])) + .ok_or_else(|| format!("reference point `{rp}` has no node")); + } + Target::Set(s) => r.sets.get(s).unwrap_or(&empty), + }; + let name = target.set().unwrap_or_default(); + if set.nodes.is_empty() { + return Err(format!("set `{name}` has no nodes")); + } + Ok(Some(match distribution { + Distribution::PerNode => set.nodes.iter().map(|n| (*n, vector)).collect(), + Distribution::Total if set.faces.is_empty() => { + let share = 1.0 / set.nodes.len() as f64; + set.nodes.iter().map(|n| (*n, scale(vector, share))).collect() + } + Distribution::Total => { + let w = self.surface_weights(&set.faces).map_err(|e| e.to_string())?; + let area: f64 = w.values().sum(); + if !(area > 0.0 && area.is_finite()) { + return Err(format!("the faces of set `{name}` have no area to share the load by")); + } + w.into_iter().map(|(n, a)| (n, scale(vector, a / area))).collect() + } + })) + } + + /// Whether a vector is parallel to the axis of a two-node element (or the element has no + /// length). + fn parallel_to_axis(&self, nodes: &[u32], v: [f64; 3]) -> bool { + let [Some(a), Some(b)] = [nodes.first(), nodes.get(1)].map(|n| n.and_then(|n| self.nodes.get(n))) else { return false }; + let t = [b.xyz[0] - a.xyz[0], b.xyz[1] - a.xyz[1], b.xyz[2] - a.xyz[2]]; + let c = [t[1] * v[2] - t[2] * v[1], t[2] * v[0] - t[0] * v[2], t[0] * v[1] - t[1] * v[0]]; + let norm = |x: [f64; 3]| (x[0] * x[0] + x[1] * x[1] + x[2] * x[2]).sqrt(); + !(norm(c) > 1e-9 * norm(t) * norm(v)) + } +} diff --git a/crates/model/src/constraints.rs b/crates/model/src/constraints.rs new file mode 100644 index 0000000..764d2b2 --- /dev/null +++ b/crates/model/src/constraints.rs @@ -0,0 +1,145 @@ +//! Reference points and couplings. +//! +//! A reference point is a single node off the mesh with three translations and three rotations. +//! A coupling ties it to a set's nodes, so a force, moment or boundary condition on the point acts +//! on a whole face (a bolt hole, a mounting pad) without hand-picked nodes. +//! +//! Each reference point keeps the node label it is written under in every deck. Unless one is +//! given, it gets the highest free label at the top of the label range when it is created, far +//! above the labels meshers use (they count up from 1), so remeshing, adding nodes and adding or +//! removing other reference points never renumber it. A mesh node on the same label is an error +//! of the check. + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::bcs::Dofs; +use crate::{Error, MAX_LABEL, Model, Result, check_coordinates, context, insert_new, not_in_use, remove_existing, replace_existing, require}; + +/// A reference point. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct ReferencePoint { + pub xyz: [f64; 3], + /// Its node label in decks; 0 when creating means "choose one" (see the module docs). + #[serde(default)] + pub label: u32, +} + +/// How a coupling ties the reference point to the set's nodes. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum CouplingKind { + /// The nodes follow the reference point rigidly in the coupled degrees of freedom. + Kinematic, + /// The reference point's load is spread over the nodes (weighted equally), and its motion is + /// their weighted average; the region stays free to deform. + Distributing, +} + +/// A coupling between a reference point and a set's nodes. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Coupling { + pub reference_point: String, + pub set: String, + pub kind: CouplingKind, + /// The coupled degrees of freedom (all six by default). + #[serde(default = "Dofs::all")] + pub dofs: Dofs, +} + +impl Model { + pub fn reference_points(&self) -> &BTreeMap { + &self.reference_points + } + + pub fn reference_point(&self, name: &str) -> Option<&ReferencePoint> { + self.reference_points.get(name) + } + + /// Check a reference point's coordinates and label, choosing a label when it has none. + fn place_reference_point(&self, name: &str, mut point: ReferencePoint, current: Option) -> Result { + check_coordinates(&point.xyz).map_err(|e| context("reference point", name, e))?; + let taken = |l: u32| self.reference_points.iter().any(|(n, p)| n != name && p.label == l); + if point.label == 0 { + point.label = match current { + Some(l) => l, + None => (1..=MAX_LABEL) + .rev() + .take(crate::MAX_DEFINITIONS + 1) + .find(|l| !taken(*l)) + .ok_or_else(|| Error::Limit(format!("no node label left for reference point `{name}`")))?, + }; + } + if point.label > MAX_LABEL { + return Err(context("reference point", name, Error::Invalid(format!("label {} is outside 1..={MAX_LABEL}", point.label)))); + } + if let Some((other, _)) = self.reference_points.iter().find(|(n, p)| *n != name && p.label == point.label) { + return Err(context("reference point", name, Error::Invalid(format!("label {} is taken by reference point `{other}`", point.label)))); + } + Ok(point) + } + + pub fn add_reference_point(&mut self, name: &str, point: ReferencePoint) -> Result<()> { + let point = self.place_reference_point(name, point, None)?; + insert_new(&mut self.reference_points, "reference point", name, point) + } + + /// Replace a reference point; a label of 0 keeps its current one. + pub fn replace_reference_point(&mut self, name: &str, point: ReferencePoint) -> Result<()> { + let current = self.reference_points.get(name).map(|p| p.label).ok_or_else(|| Error::NotFound(format!("reference point `{name}`")))?; + let point = self.place_reference_point(name, point, Some(current))?; + replace_existing(&mut self.reference_points, "reference point", name, point) + } + + pub fn remove_reference_point(&mut self, name: &str) -> Result<()> { + let mut users: Vec = + self.couplings.iter().filter(|(_, c)| c.reference_point == name).map(|(n, _)| format!("coupling `{n}`")).collect(); + users.extend( + self.boundary_conditions + .iter() + .filter(|(_, b)| b.target().reference_point() == Some(name)) + .map(|(n, _)| format!("boundary condition `{n}`")), + ); + users.extend(self.loads.iter().filter(|(_, l)| l.reference_point() == Some(name)).map(|(n, _)| format!("load `{n}`"))); + not_in_use("reference point", name, users)?; + remove_existing(&mut self.reference_points, "reference point", name).map(|_| ()) + } + + /// The node label of each reference point in decks. + pub fn reference_point_nodes(&self) -> BTreeMap { + self.reference_points.iter().map(|(n, p)| (n.clone(), p.label)).collect() + } + + pub fn couplings(&self) -> &BTreeMap { + &self.couplings + } + + pub fn coupling(&self, name: &str) -> Option<&Coupling> { + self.couplings.get(name) + } + + fn check_coupling(&self, name: &str, c: &Coupling) -> Result<()> { + if c.dofs.is_empty() { + return Err(context("coupling", name, Error::Invalid("couples no degree of freedom".into()))); + } + require(&self.reference_points, "reference point", &c.reference_point).map_err(|e| context("coupling", name, e))?; + require(&self.sets, "set", &c.set).map_err(|e| context("coupling", name, e)) + } + + pub fn add_coupling(&mut self, name: &str, coupling: Coupling) -> Result<()> { + self.check_coupling(name, &coupling)?; + insert_new(&mut self.couplings, "coupling", name, coupling) + } + + pub fn replace_coupling(&mut self, name: &str, coupling: Coupling) -> Result<()> { + self.check_coupling(name, &coupling)?; + replace_existing(&mut self.couplings, "coupling", name, coupling) + } + + pub fn remove_coupling(&mut self, name: &str) -> Result<()> { + remove_existing(&mut self.couplings, "coupling", name).map(|_| ()) + } +} diff --git a/crates/model/src/geometry.rs b/crates/model/src/geometry.rs new file mode 100644 index 0000000..f567736 --- /dev/null +++ b/crates/model/src/geometry.rs @@ -0,0 +1,360 @@ +//! How the mesh lies on the geometry. +//! +//! Geometry is referred to by name only ([`GeometryRef`]): a body's name and, inside it, the +//! persistent name of a face, edge or vertex (for SolidCraft designs, names such as `F2:side:l4` +//! that follow the face through edits of the design). The model never sees the geometry itself. +//! +//! A mesher records where it put each mesh entity ([`Classification`]): every element on the +//! body (solids), face (shells) or edge (beams) it fills, every boundary element face on the +//! geometry face it lies on, element edges on geometry edges and nodes on geometry vertices. +//! Sets defined by geometry names resolve through these tags ([`crate::sets`]), so they follow a +//! remesh: the new mesh carries new tags for the same names. + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::{Error, Model, Result}; + +/// Longest body or entity name accepted (persistent names can be long; they are not written to +/// decks as such). +pub const MAX_GEOMETRY_NAME_LEN: usize = 256; + +/// Most distinct geometry entities a mesh can be tagged with. +pub const MAX_ENTITIES: usize = 1_000_000; + +/// An entity of a body. +#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub enum Entity { + /// The whole body. + Body, + Face(String), + Edge(String), + Vertex(String), +} + +impl Entity { + /// `body`, `face`, `edge` or `vertex`. + pub fn kind(&self) -> &'static str { + match self { + Entity::Body => "body", + Entity::Face(_) => "face", + Entity::Edge(_) => "edge", + Entity::Vertex(_) => "vertex", + } + } + + /// The entity's name within its body (`None` for the body itself). + pub fn name(&self) -> Option<&str> { + match self { + Entity::Body => None, + Entity::Face(n) | Entity::Edge(n) | Entity::Vertex(n) => Some(n), + } + } + + /// The topological dimension: 3 body, 2 face, 1 edge, 0 vertex. + pub fn dimension(&self) -> u8 { + match self { + Entity::Body => 3, + Entity::Face(_) => 2, + Entity::Edge(_) => 1, + Entity::Vertex(_) => 0, + } + } +} + +/// A named piece of geometry: `{"body": "Body1"}`, `{"body": "Body1", "face": "F2:top"}`, or +/// with `edge` or `vertex` in place of `face`. +#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(try_from = "GeometryRefRepr", into = "GeometryRefRepr")] +pub struct GeometryRef { + pub body: String, + pub entity: Entity, +} + +#[derive(Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +struct GeometryRefRepr { + body: String, + #[serde(default, skip_serializing_if = "Option::is_none")] + face: Option, + #[serde(default, skip_serializing_if = "Option::is_none")] + edge: Option, + #[serde(default, skip_serializing_if = "Option::is_none")] + vertex: Option, +} + +impl TryFrom for GeometryRef { + type Error = Error; + fn try_from(r: GeometryRefRepr) -> Result { + let entity = match (r.face, r.edge, r.vertex) { + (None, None, None) => Entity::Body, + (Some(n), None, None) => Entity::Face(n), + (None, Some(n), None) => Entity::Edge(n), + (None, None, Some(n)) => Entity::Vertex(n), + _ => return Err(Error::Invalid("a geometry reference names at most one of `face`, `edge` or `vertex`".into())), + }; + GeometryRef::new(r.body, entity) + } +} + +impl From for GeometryRefRepr { + fn from(g: GeometryRef) -> Self { + let mut r = GeometryRefRepr { body: g.body, face: None, edge: None, vertex: None }; + match g.entity { + Entity::Body => {} + Entity::Face(n) => r.face = Some(n), + Entity::Edge(n) => r.edge = Some(n), + Entity::Vertex(n) => r.vertex = Some(n), + } + r + } +} + +fn check_geometry_name(what: &str, name: &str) -> Result<()> { + if name.is_empty() || name.len() > MAX_GEOMETRY_NAME_LEN || name.chars().any(char::is_control) { + let shown: String = name.chars().take(64).filter(|c| !c.is_control()).collect(); + return Err(Error::Invalid(format!("{what} name `{shown}` must be 1 to {MAX_GEOMETRY_NAME_LEN} characters with no control characters"))); + } + Ok(()) +} + +impl GeometryRef { + pub fn new(body: impl Into, entity: Entity) -> Result { + let body = body.into(); + check_geometry_name("body", &body)?; + if let Some(n) = entity.name() { + check_geometry_name(entity.kind(), n)?; + } + Ok(GeometryRef { body, entity }) + } + + pub fn body(body: impl Into) -> Result { + GeometryRef::new(body, Entity::Body) + } + + pub fn face(body: impl Into, face: impl Into) -> Result { + GeometryRef::new(body, Entity::Face(face.into())) + } + + /// Whether this reference (from a set) selects the tagged entity `tag`: the same body and + /// kind, and the same name or the name of one of its pieces. A B-rep face the kernel splits + /// into pieces (a full cylinder is two half-cylinders) carries the face's name with `#0`, + /// `#1`, … appended; a reference to the plain name selects all of them. + pub fn selects(&self, tag: &GeometryRef) -> bool { + if self.body != tag.body || self.entity.kind() != tag.entity.kind() { + return false; + } + match (self.entity.name(), tag.entity.name()) { + (Some(want), Some(have)) => { + have == want + || have + .strip_prefix(want) + .and_then(|r| r.strip_prefix('#')) + .is_some_and(|k| !k.is_empty() && k.bytes().all(|b| b.is_ascii_digit())) + } + _ => true, + } + } +} + +impl std::fmt::Display for GeometryRef { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self.entity.name() { + None => write!(f, "body {}", self.body), + Some(n) => write!(f, "{} {n} of body {}", self.entity.kind(), self.body), + } + } +} + +/// The geometry tags of one element. +#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)] +pub struct ElementTags { + /// The entity the element fills: a body for solids, a face for shells, an edge for beams. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub region: Option, + /// Face number (1-based, [`crate::ElementKind::faces`]) → the geometry face it lies on. + #[serde(default, skip_serializing_if = "BTreeMap::is_empty")] + pub faces: BTreeMap, + /// Edge number (1-based, [`crate::ElementKind::edges`]) → the geometry edge it lies on. + #[serde(default, skip_serializing_if = "BTreeMap::is_empty")] + pub edges: BTreeMap, +} + +/// Where the mesh lies on the geometry. Tags are indices into [`Classification::entities`]. +/// +/// It is read back from a file through [`Model`]'s deserializer, which replays every tag through +/// the checked `tag_*` methods. +#[derive(Debug, Clone, Default, PartialEq, Serialize)] +pub struct Classification { + entities: Vec, + /// Entity → tag, so interning is a lookup. + #[serde(skip)] + index: BTreeMap, + elements: BTreeMap, + /// Node → the geometry vertex it lies on. + nodes: BTreeMap, +} + +/// A [`Classification`] as stored, before its tags are checked against the mesh. +#[derive(Debug, Default, Deserialize)] +#[serde(deny_unknown_fields)] +pub(crate) struct ClassificationData { + #[serde(default)] + pub entities: Vec, + #[serde(default)] + pub elements: BTreeMap, + #[serde(default)] + pub nodes: BTreeMap, +} + +impl Classification { + /// The tagged entities; a tag is an index into this list. + pub fn entities(&self) -> &[GeometryRef] { + &self.entities + } + + pub fn entity(&self, tag: u32) -> Option<&GeometryRef> { + self.entities.get(tag as usize) + } + + pub fn elements(&self) -> &BTreeMap { + &self.elements + } + + pub fn nodes(&self) -> &BTreeMap { + &self.nodes + } + + pub fn is_empty(&self) -> bool { + self.elements.is_empty() && self.nodes.is_empty() + } + + /// The tags of every entity `reference` selects ([`GeometryRef::selects`]). + pub fn tags_selected_by(&self, reference: &GeometryRef) -> Vec { + self.entities.iter().enumerate().filter(|(_, e)| reference.selects(e)).filter_map(|(i, _)| u32::try_from(i).ok()).collect() + } + + pub(crate) fn forget_element(&mut self, label: u32) { + self.elements.remove(&label); + } + + pub(crate) fn forget_node(&mut self, label: u32) { + self.nodes.remove(&label); + } +} + +impl Model { + pub fn classification(&self) -> &Classification { + &self.classification + } + + /// The tag for a geometry entity, adding it on first use. Meshers call this once per entity, + /// then tag mesh entities with the number. + pub fn geometry_tag(&mut self, entity: &GeometryRef) -> Result { + let c = &mut self.classification; + if let Some(&tag) = c.index.get(entity) { + return Ok(tag); + } + if c.entities.len() >= MAX_ENTITIES { + return Err(Error::Limit(format!("a mesh is tagged with at most {MAX_ENTITIES} geometry entities"))); + } + let tag = u32::try_from(c.entities.len()).map_err(|_| Error::Limit("too many geometry entities".into()))?; + c.entities.push(entity.clone()); + c.index.insert(entity.clone(), tag); + Ok(tag) + } + + /// Check a tag names an entity of the given dimension. + fn check_tag(&self, tag: u32, dimension: u8) -> Result<()> { + let e = self.classification.entity(tag).ok_or_else(|| Error::NotFound(format!("geometry tag {tag}")))?; + if e.entity.dimension() != dimension { + return Err(Error::Invalid(format!("{e} is a {}, not a {}-dimensional entity", e.entity.kind(), dimension))); + } + Ok(()) + } + + fn element_kind(&self, element: u32) -> Result { + Ok(self.element(element).ok_or_else(|| Error::NotFound(format!("element {element}")))?.kind) + } + + /// Check that `number` (1-based) is one of `parts` of an element, naming `what` in the error. + fn check_part(element: u32, kind: crate::ElementKind, what: &str, number: u8, parts: usize) -> Result<()> { + if number == 0 || usize::from(number) > parts { + return Err(Error::Invalid(format!("element {element} ({}) has no {what} {number}", kind.name()))); + } + Ok(()) + } + + /// Record that an element fills a body (solids), a face (shells) or an edge (line elements). + pub fn tag_element(&mut self, element: u32, tag: u32) -> Result<()> { + let kind = self.element_kind(element)?; + self.check_tag(tag, kind.dimension())?; + self.classification.elements.entry(element).or_default().region = Some(tag); + Ok(()) + } + + /// Record that face `face` (1-based) of an element lies on a geometry face. + pub fn tag_element_face(&mut self, element: u32, face: u8, tag: u32) -> Result<()> { + let kind = self.element_kind(element)?; + Self::check_part(element, kind, "face", face, kind.faces().len())?; + self.check_tag(tag, 2)?; + self.classification.elements.entry(element).or_default().faces.insert(face, tag); + Ok(()) + } + + /// Record that edge `edge` (1-based) of an element lies on a geometry edge. + pub fn tag_element_edge(&mut self, element: u32, edge: u8, tag: u32) -> Result<()> { + let kind = self.element_kind(element)?; + Self::check_part(element, kind, "edge", edge, kind.edges().len())?; + self.check_tag(tag, 1)?; + self.classification.elements.entry(element).or_default().edges.insert(edge, tag); + Ok(()) + } + + /// Record that a node lies on a geometry vertex. + pub fn tag_node(&mut self, node: u32, tag: u32) -> Result<()> { + if self.node(node).is_none() { + return Err(Error::NotFound(format!("node {node}"))); + } + self.check_tag(tag, 0)?; + self.classification.nodes.insert(node, tag); + Ok(()) + } + + /// Drop every geometry tag (before a remesh writes new ones). + pub fn clear_geometry_tags(&mut self) { + self.classification = Classification::default(); + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn references_parse_and_select_split_faces() { + let r: GeometryRef = serde_json::from_value(serde_json::json!({"body": "Body1", "face": "F2:side"})).unwrap(); + assert_eq!(r, GeometryRef::face("Body1", "F2:side").unwrap()); + assert_eq!(serde_json::to_value(&r).unwrap(), serde_json::json!({"body": "Body1", "face": "F2:side"})); + assert!(r.selects(&GeometryRef::face("Body1", "F2:side").unwrap())); + assert!(r.selects(&GeometryRef::face("Body1", "F2:side#1").unwrap())); + assert!(!r.selects(&GeometryRef::face("Body1", "F2:side#").unwrap())); + assert!(!r.selects(&GeometryRef::face("Body1", "F2:sides").unwrap())); + assert!(!r.selects(&GeometryRef::face("Body2", "F2:side").unwrap())); + assert!(!r.selects(&GeometryRef::new("Body1", Entity::Edge("F2:side".into())).unwrap())); + assert!(GeometryRef::body("Body1").unwrap().selects(&GeometryRef::body("Body1").unwrap())); + + for bad in [ + serde_json::json!({"body": ""}), + serde_json::json!({"body": "B", "face": "a", "edge": "b"}), + serde_json::json!({"body": "B", "face": ""}), + serde_json::json!({"body": "B", "face": "a\nb"}), + serde_json::json!({"body": "x".repeat(MAX_GEOMETRY_NAME_LEN + 1)}), + serde_json::json!({"body": "B", "colour": "red"}), + serde_json::json!({"face": "a"}), + ] { + assert!(serde_json::from_value::(bad.clone()).is_err(), "{bad}"); + } + } +} diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index e4454b2..da05b00 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -1,23 +1,64 @@ //! The solver-neutral finite element model. //! -//! A [`Model`] holds what every solver deck describes: nodes and elements addressed by the -//! positive integer labels the decks use. Materials, sections, sets, boundary conditions, loads -//! and steps join it as their tasks land (plan/architecture.md §4). Solver formats (Abaqus input, -//! Nastran bulk data, LS-DYNA keyword) translate to and from this model; nothing here knows any -//! of them. +//! A [`Model`] holds what every solver deck describes, in no solver's vocabulary: +//! - the mesh: nodes and elements addressed by the positive integer labels the decks use +//! ([`mesh`]), and how the mesh lies on the geometry ([`geometry`]: each element, element face, +//! element edge and node tagged with the persistent name of the body, face, edge or vertex it +//! was meshed on); +//! - named sets ([`sets`]) defined by geometry names or explicit labels, resolved to nodes, +//! elements and element faces against the current mesh, so they survive remeshing; +//! - materials and the unit system ([`materials`], [`units`]), sections with their element +//! formulation ([`sections`]), reference points and couplings ([`constraints`]); +//! - boundary conditions ([`bcs`]), loads ([`loads`]) and steps with an explicit table of which +//! boundary conditions and loads are active in each step ([`steps`]); +//! - the check that names every missing or bad item, and the resolved view the deck writers and +//! the solver read ([`check`]). +//! +//! Solver formats (Abaqus input, Nastran bulk data, LS-DYNA keyword) translate to and from this +//! model; nothing here knows any of them. //! //! The model is unit-consistent: it stores numbers in whatever consistent system the user chose -//! (for example mm, tonne, s, N, MPa) and never converts them. +//! (for example mm, tonne, s, N, MPa), never converts them, and only warns when a magnitude looks +//! wrong for the chosen system. //! -//! Every mutation validates its input, because labels and coordinates come from files, scripts -//! and the control channel: a bad value is an [`Error`], never a panic. +//! Every mutation validates its input, because labels, names and numbers come from files, scripts +//! and the control channel: a bad value is an [`Error`], never a panic. Definitions that refer to +//! other definitions (a section to its set and material, a load to its set) must refer to ones +//! that exist, and a definition still referred to cannot be deleted. What depends on the mesh (a +//! set's geometry names and labels) is checked when the model is resolved, because remeshing +//! replaces the mesh under the definitions. #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] #![forbid(unsafe_code)] +pub mod bcs; +pub mod check; +pub mod constraints; +pub mod geometry; +pub mod loads; +pub mod materials; +mod mesh; +pub mod nodal; +pub mod sections; +pub mod sets; +pub mod steps; +pub mod units; + use std::collections::BTreeMap; use serde::{Deserialize, Serialize}; +pub use bcs::{BoundaryCondition, Dofs, Symmetry, Target}; +pub use check::{Diagnostic, Resolved, ResolvedSet, Severity}; +pub use constraints::{Coupling, CouplingKind, ReferencePoint}; +pub use geometry::{Classification, Entity, GeometryRef}; +pub use loads::{Distribution, Load}; +pub use materials::{Elastic, Material}; +pub use mesh::{Element, ElementFace, ElementKind, Node, check_coordinates}; +pub use sections::{BeamProfile, Formulation, Integration, Section, SectionKind}; +pub use sets::Set; +pub use steps::{OutputVariable, Procedure, Step}; +pub use units::UnitSystem; + /// The largest node or element label. Nastran's GRID and element ids stop at 99 999 999 (eight /// digits in a small field), the tightest of the three formats. pub const MAX_LABEL: u32 = 99_999_999; @@ -29,6 +70,13 @@ pub const MAX_ELEMENTS: usize = 50_000_000; /// Coordinates beyond this magnitude are rejected (1e9 model units). pub const MAX_COORD: f64 = 1e9; +/// Most definitions of one kind (sets, materials, sections, …) a model holds. +pub const MAX_DEFINITIONS: usize = 100_000; + +/// Longest name of a definition. Names are ASCII: a letter, then letters, digits, `_` or `-`, so +/// every deck format can carry them unquoted. +pub const MAX_NAME_LEN: usize = 64; + #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum Error { #[error("invalid input: {0}")] @@ -37,267 +85,214 @@ pub enum Error { NotFound(String), #[error("limit reached: {0}")] Limit(String), + /// A definition cannot be deleted while another refers to it. + #[error("in use: {0}")] + InUse(String), } pub type Result = std::result::Result; -/// A node: a point in model space. -#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] -pub struct Node { - pub xyz: [f64; 3], -} - -/// Element topology. Node order is corners first, then mid-side nodes, with positive volume for -/// a right-handed ordering: -/// - tet10: mid-side nodes 5–10 on edges 1-2, 2-3, 3-1, 1-4, 2-4, 3-4; -/// - wedge6: triangle 1-2-3, then 4-5-6 on the side its right-hand normal points to; -/// - hex8/hex20: face 1-2-3-4, then 5-6-7-8 above it; hex20's mid-side nodes 9–12 lie on the -/// edges of 1-2-3-4, 13–16 on those of 5-6-7-8, and 17–20 on edges 1-5, 2-6, 3-7, 4-8. +/// The finite element model. /// -/// This is the Abaqus (and VTK) order. Formats that order nodes differently, such as Nastran's -/// CHEXA with 20 nodes (vertical edges before the top face's edges), are reordered by their deck -/// readers and writers; the model never stores another order. -#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] -#[serde(rename_all = "snake_case")] -pub enum ElementKind { - /// Two-node axial member. - Rod2, - /// Two-node beam. - Beam2, - /// Three-node triangle (shell or membrane). - Tri3, - /// Six-node triangle. - Tri6, - /// Four-node quadrilateral. - Quad4, - /// Eight-node quadrilateral. - Quad8, - /// Four-node tetrahedron. - Tet4, - /// Ten-node tetrahedron. - Tet10, - /// Six-node wedge (pentahedron). - Wedge6, - /// Eight-node hexahedron. - Hex8, - /// Twenty-node hexahedron. - Hex20, -} - -impl ElementKind { - pub const ALL: [ElementKind; 11] = [ - ElementKind::Rod2, - ElementKind::Beam2, - ElementKind::Tri3, - ElementKind::Tri6, - ElementKind::Quad4, - ElementKind::Quad8, - ElementKind::Tet4, - ElementKind::Tet10, - ElementKind::Wedge6, - ElementKind::Hex8, - ElementKind::Hex20, - ]; - - pub fn node_count(self) -> usize { - match self { - ElementKind::Rod2 | ElementKind::Beam2 => 2, - ElementKind::Tri3 => 3, - ElementKind::Quad4 | ElementKind::Tet4 => 4, - ElementKind::Tri6 | ElementKind::Wedge6 => 6, - ElementKind::Quad8 | ElementKind::Hex8 => 8, - ElementKind::Tet10 => 10, - ElementKind::Hex20 => 20, - } - } - - /// Topological dimension: 1 for line elements, 2 for shells, 3 for solids. - pub fn dimension(self) -> u8 { - match self { - ElementKind::Rod2 | ElementKind::Beam2 => 1, - ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => 2, - ElementKind::Tet4 | ElementKind::Tet10 | ElementKind::Wedge6 | ElementKind::Hex8 | ElementKind::Hex20 => 3, - } - } - - /// The name used in commands and `.modelcraft` files (`hex8`, `tet10`, …). - pub fn name(self) -> &'static str { - match self { - ElementKind::Rod2 => "rod2", - ElementKind::Beam2 => "beam2", - ElementKind::Tri3 => "tri3", - ElementKind::Tri6 => "tri6", - ElementKind::Quad4 => "quad4", - ElementKind::Quad8 => "quad8", - ElementKind::Tet4 => "tet4", - ElementKind::Tet10 => "tet10", - ElementKind::Wedge6 => "wedge6", - ElementKind::Hex8 => "hex8", - ElementKind::Hex20 => "hex20", - } - } - - pub fn from_name(name: &str) -> Option { - ElementKind::ALL.into_iter().find(|k| k.name() == name) - } -} - -/// An element: a topology and its node labels. -#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] -pub struct Element { - pub kind: ElementKind, - pub nodes: Vec, -} - -/// Nodes and elements by label. +/// It deserializes by replaying what it reads through the same checked methods commands use, so +/// a model read from a file holds the same invariants as one built by commands. #[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)] +#[serde(try_from = "ModelData")] pub struct Model { nodes: BTreeMap, elements: BTreeMap, + classification: Classification, + units: UnitSystem, + sets: BTreeMap, + materials: BTreeMap, + sections: BTreeMap, + reference_points: BTreeMap, + couplings: BTreeMap, + boundary_conditions: BTreeMap, + loads: BTreeMap, + /// In analysis order. + steps: Vec, } -/// Check that every coordinate is a finite number within ±[`MAX_COORD`] (the rule for nodes, -/// shared by everything that takes coordinates from outside). -pub fn check_coordinates(xyz: &[f64]) -> Result<()> { - match xyz.iter().find(|c| !(c.abs() <= MAX_COORD)) { - Some(c) => Err(Error::Invalid(format!("coordinate {c} is not a finite number within ±{MAX_COORD}"))), - None => Ok(()), +/// A [`Model`] as stored, before it is checked. +#[derive(Deserialize)] +#[serde(deny_unknown_fields)] +struct ModelData { + #[serde(default)] + nodes: BTreeMap, + #[serde(default)] + elements: BTreeMap, + #[serde(default)] + classification: geometry::ClassificationData, + #[serde(default)] + units: UnitSystem, + #[serde(default)] + sets: BTreeMap, + #[serde(default)] + materials: BTreeMap, + #[serde(default)] + sections: BTreeMap, + #[serde(default)] + reference_points: BTreeMap, + #[serde(default)] + couplings: BTreeMap, + #[serde(default)] + boundary_conditions: BTreeMap, + #[serde(default)] + loads: BTreeMap, + #[serde(default)] + steps: Vec, +} + +impl TryFrom for Model { + type Error = Error; + + fn try_from(d: ModelData) -> Result { + let mut m = Model { units: d.units, ..Model::default() }; + for (label, node) in d.nodes { + m.set_node(label, node.xyz)?; + } + for (label, e) in d.elements { + m.set_element(label, e.kind, e.nodes)?; + } + for (i, entity) in d.classification.entities.iter().enumerate() { + if usize::try_from(m.geometry_tag(entity)?).ok() != Some(i) { + return Err(Error::Invalid(format!("{entity} is listed twice among the geometry tags"))); + } + } + for (element, tags) in d.classification.elements { + if let Some(t) = tags.region { + m.tag_element(element, t)?; + } + for (face, t) in tags.faces { + m.tag_element_face(element, face, t)?; + } + for (edge, t) in tags.edges { + m.tag_element_edge(element, edge, t)?; + } + } + for (node, t) in d.classification.nodes { + m.tag_node(node, t)?; + } + for (n, v) in d.sets { + m.add_set(&n, v)?; + } + for (n, v) in d.materials { + m.add_material(&n, v)?; + } + for (n, v) in d.sections { + m.add_section(&n, v)?; + } + for (n, v) in d.reference_points { + m.add_reference_point(&n, v)?; + } + for (n, v) in d.couplings { + m.add_coupling(&n, v)?; + } + for (n, v) in d.boundary_conditions { + m.add_boundary_condition(&n, v)?; + } + for (n, v) in d.loads { + m.add_load(&n, v)?; + } + for s in d.steps { + m.add_step(s, None)?; + } + Ok(m) } } -fn check_label(what: &str, label: u32) -> Result<()> { - if label == 0 || label > MAX_LABEL { - return Err(Error::Invalid(format!("{what} label {label} is outside 1..={MAX_LABEL}"))); +/// Check a definition's name: 1..=[`MAX_NAME_LEN`] ASCII characters, a letter first, then +/// letters, digits, `_` or `-`. +pub fn check_name(what: &str, name: &str) -> Result<()> { + let mut chars = name.chars(); + let ok = name.len() <= MAX_NAME_LEN + && chars.next().is_some_and(|c| c.is_ascii_alphabetic()) + && chars.all(|c| c.is_ascii_alphanumeric() || c == '_' || c == '-'); + if ok { + Ok(()) + } else { + let shown: String = name.chars().take(MAX_NAME_LEN).collect(); + Err(Error::Invalid(format!("{what} name `{shown}` must be 1 to {MAX_NAME_LEN} characters: a letter, then letters, digits, `_` or `-`"))) } +} + +/// Check that a value is a finite number, naming it in the error. +pub(crate) fn finite(what: &str, v: f64) -> Result { + if v.is_finite() { Ok(v) } else { Err(Error::Invalid(format!("{what} must be a finite number, got {v}"))) } +} + +/// Check that a value is a finite number above zero, naming it in the error. +pub(crate) fn positive(what: &str, v: f64) -> Result { + if v > 0.0 && v.is_finite() { Ok(v) } else { Err(Error::Invalid(format!("{what} must be a positive number, got {v}"))) } +} + +/// Add a new definition under a fresh, valid name. +fn insert_new(map: &mut BTreeMap, what: &str, name: &str, value: T) -> Result<()> { + check_name(what, name)?; + if map.contains_key(name) { + return Err(Error::Invalid(format!("{what} `{name}` already exists"))); + } + if map.len() >= MAX_DEFINITIONS { + return Err(Error::Limit(format!("a model holds at most {MAX_DEFINITIONS} of each kind of definition"))); + } + map.insert(name.to_string(), value); Ok(()) } +/// Replace an existing definition. +fn replace_existing(map: &mut BTreeMap, what: &str, name: &str, value: T) -> Result<()> { + let slot = map.get_mut(name).ok_or_else(|| Error::NotFound(format!("{what} `{name}`")))?; + *slot = value; + Ok(()) +} + +/// Prefix an error with the definition it is about. +pub(crate) fn context(what: &str, name: &str, e: Error) -> Error { + let shown: String = name.chars().take(MAX_NAME_LEN).collect(); + match e { + Error::Invalid(m) => Error::Invalid(format!("{what} `{shown}`: {m}")), + Error::NotFound(m) => Error::NotFound(format!("{what} `{shown}`: {m}")), + Error::Limit(m) => Error::Limit(format!("{what} `{shown}`: {m}")), + Error::InUse(m) => Error::InUse(format!("{what} `{shown}`: {m}")), + } +} + +/// Remove an existing definition. +fn remove_existing(map: &mut BTreeMap, what: &str, name: &str) -> Result { + map.remove(name).ok_or_else(|| Error::NotFound(format!("{what} `{name}`"))) +} + +/// Fail with [`Error::InUse`] when any referrer is listed. +fn not_in_use(what: &str, name: &str, referrers: Vec) -> Result<()> { + if referrers.is_empty() { Ok(()) } else { Err(Error::InUse(format!("{what} `{name}` is used by {}", referrers.join(", ")))) } +} + +fn require(map: &BTreeMap, what: &str, name: &str) -> Result<()> { + if map.contains_key(name) { Ok(()) } else { Err(Error::NotFound(format!("{what} `{name}`"))) } +} + impl Model { - pub fn nodes(&self) -> &BTreeMap { - &self.nodes + pub fn units(&self) -> UnitSystem { + self.units } - pub fn elements(&self) -> &BTreeMap { - &self.elements + /// Label the unit system the model's numbers are in (nothing is converted). + pub fn set_units(&mut self, units: UnitSystem) { + self.units = units; } - pub fn node(&self, label: u32) -> Option<&Node> { - self.nodes.get(&label) - } - - pub fn element(&self, label: u32) -> Option<&Element> { - self.elements.get(&label) - } - - /// Add a node, or move an existing one (decks redefine nodes the same way). - pub fn set_node(&mut self, label: u32, xyz: [f64; 3]) -> Result<()> { - check_label("node", label)?; - check_coordinates(&xyz).map_err(|e| match e { - Error::Invalid(m) => Error::Invalid(format!("node {label}: {m}")), - e => e, - })?; - if !self.nodes.contains_key(&label) && self.nodes.len() >= MAX_NODES { - return Err(Error::Limit(format!("a model holds at most {MAX_NODES} nodes"))); - } - self.nodes.insert(label, Node { xyz }); - Ok(()) - } - - /// Add an element, or replace an existing one. Its nodes must exist and be distinct. - pub fn set_element(&mut self, label: u32, kind: ElementKind, nodes: Vec) -> Result<()> { - check_label("element", label)?; - if nodes.len() != kind.node_count() { - return Err(Error::Invalid(format!("element {label}: {} needs {} nodes, got {}", kind.name(), kind.node_count(), nodes.len()))); - } - for (i, n) in nodes.iter().enumerate() { - if !self.nodes.contains_key(n) { - return Err(Error::NotFound(format!("element {label}: node {n}"))); - } - if nodes.get(..i).is_some_and(|before| before.contains(n)) { - return Err(Error::Invalid(format!("element {label}: node {n} appears twice"))); - } - } - if !self.elements.contains_key(&label) && self.elements.len() >= MAX_ELEMENTS { - return Err(Error::Limit(format!("a model holds at most {MAX_ELEMENTS} elements"))); - } - self.elements.insert(label, Element { kind, nodes }); - Ok(()) - } - - /// Remove a node that no element uses. - pub fn remove_node(&mut self, label: u32) -> Result<()> { - if let Some((e, _)) = self.elements.iter().find(|(_, e)| e.nodes.contains(&label)) { - return Err(Error::Invalid(format!("node {label} is used by element {e}"))); - } - self.nodes.remove(&label).map(|_| ()).ok_or_else(|| Error::NotFound(format!("node {label}"))) - } - - pub fn remove_element(&mut self, label: u32) -> Result<()> { - self.elements.remove(&label).map(|_| ()).ok_or_else(|| Error::NotFound(format!("element {label}"))) - } - - /// Element counts by kind, in [`ElementKind::ALL`] order, kinds with none left out. - pub fn element_counts(&self) -> Vec<(ElementKind, usize)> { - ElementKind::ALL.into_iter().map(|k| (k, self.elements.values().filter(|e| e.kind == k).count())).filter(|(_, n)| *n > 0).collect() - } -} - -#[cfg(test)] -mod tests { - use super::*; - - fn cube() -> Model { - let mut m = Model::default(); - let corners = [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]]; - for (i, c) in corners.into_iter().enumerate() { - m.set_node(i as u32 + 1, c).unwrap(); - } - m.set_element(1, ElementKind::Hex8, (1..=8).collect()).unwrap(); - m - } - - #[test] - fn builds_a_hex() { - let m = cube(); - assert_eq!(m.nodes().len(), 8); - assert_eq!(m.element_counts(), vec![(ElementKind::Hex8, 1)]); - } - - #[test] - fn rejects_bad_labels_and_coordinates() { - let mut m = Model::default(); - assert!(m.set_node(0, [0.; 3]).is_err()); - assert!(m.set_node(MAX_LABEL + 1, [0.; 3]).is_err()); - assert!(m.set_node(1, [f64::NAN, 0., 0.]).is_err()); - assert!(m.set_node(1, [f64::INFINITY, 0., 0.]).is_err()); - assert!(m.set_node(1, [2e9, 0., 0.]).is_err()); - assert!(m.nodes().is_empty()); - } - - #[test] - fn rejects_bad_connectivity() { - let mut m = cube(); - assert!(matches!(m.set_element(2, ElementKind::Tet4, vec![1, 2, 3]), Err(Error::Invalid(_)))); - assert!(matches!(m.set_element(2, ElementKind::Tet4, vec![1, 2, 3, 99]), Err(Error::NotFound(_)))); - assert!(matches!(m.set_element(2, ElementKind::Tet4, vec![1, 2, 3, 3]), Err(Error::Invalid(_)))); - assert_eq!(m.elements().len(), 1); - } - - #[test] - fn used_nodes_cannot_be_removed() { - let mut m = cube(); - assert!(m.remove_node(1).is_err()); - m.remove_element(1).unwrap(); - m.remove_node(1).unwrap(); - assert!(m.remove_node(1).is_err()); - } - - #[test] - fn kind_names_round_trip() { - for k in ElementKind::ALL { - assert_eq!(ElementKind::from_name(k.name()), Some(k)); - assert_eq!(serde_json::to_value(k).unwrap(), serde_json::json!(k.name())); - } - assert_eq!(ElementKind::from_name("C3D8"), None); + /// Counts of each kind of definition, for summaries. + pub fn definition_counts(&self) -> [(&'static str, usize); 8] { + [ + ("sets", self.sets.len()), + ("materials", self.materials.len()), + ("sections", self.sections.len()), + ("reference_points", self.reference_points.len()), + ("couplings", self.couplings.len()), + ("boundary_conditions", self.boundary_conditions.len()), + ("loads", self.loads.len()), + ("steps", self.steps.len()), + ] } } diff --git a/crates/model/src/loads.rs b/crates/model/src/loads.rs new file mode 100644 index 0000000..4846d55 --- /dev/null +++ b/crates/model/src/loads.rs @@ -0,0 +1,145 @@ +//! Loads. Which steps a load is active in is listed explicitly by each step ([`crate::steps`]). + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::bcs::Target; +use crate::{Error, Model, Result, context, finite, insert_new, not_in_use, remove_existing, replace_existing, require}; + +/// How a force or moment on a set is shared among the set's nodes. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum Distribution { + /// The vector is the total. On a set with element faces it is split by the faces' consistent + /// nodal weights (area-weighted, [`Model::surface_weights`]); on a set of nodes only, equally. + /// The total is the same on every mesh of the set. + #[default] + Total, + /// The vector is applied at every node of the set, so the total grows with the number of + /// nodes: it depends on the mesh. + PerNode, +} + +/// A load. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] +pub enum Load { + /// A force on the target's nodes, shared as `distribution` says (a reference point takes it + /// whole). + Force { + target: Target, + vector: [f64; 3], + #[serde(default)] + distribution: Distribution, + }, + /// A moment on the target's nodes, shared as `distribution` says. Solid elements have no + /// rotations, so a moment on a solid goes through a coupled reference point. + Moment { + target: Target, + vector: [f64; 3], + #[serde(default)] + distribution: Distribution, + }, + /// A force per area on the set's element faces, as consistent nodal forces. + Traction { set: String, vector: [f64; 3] }, + /// Uniform pressure on the set's element faces, positive pushing into the elements (against + /// the outward normal). + Pressure { set: String, magnitude: f64 }, + /// Gravity: an acceleration of the mass of the set's elements, or of the whole model when no + /// set is given. + Gravity { + acceleration: [f64; 3], + #[serde(default, skip_serializing_if = "Option::is_none")] + set: Option, + }, + /// A force per volume on the set's elements. + BodyForce { set: String, vector: [f64; 3] }, +} + +impl Load { + /// `force`, `moment`, `traction`, `pressure`, `gravity` or `body_force`. + pub fn kind(&self) -> &'static str { + match self { + Load::Force { .. } => "force", + Load::Moment { .. } => "moment", + Load::Traction { .. } => "traction", + Load::Pressure { .. } => "pressure", + Load::Gravity { .. } => "gravity", + Load::BodyForce { .. } => "body_force", + } + } + + /// The set the load acts on, if any. + pub fn set(&self) -> Option<&str> { + match self { + Load::Force { target, .. } | Load::Moment { target, .. } => target.set(), + Load::Pressure { set, .. } | Load::BodyForce { set, .. } | Load::Traction { set, .. } => Some(set), + Load::Gravity { set, .. } => set.as_deref(), + } + } + + /// The reference point the load acts on, if any. + pub fn reference_point(&self) -> Option<&str> { + match self { + Load::Force { target, .. } | Load::Moment { target, .. } => target.reference_point(), + _ => None, + } + } + + pub fn validate(&self) -> Result<()> { + let vector = match self { + Load::Force { vector, .. } | Load::Moment { vector, .. } | Load::BodyForce { vector, .. } | Load::Traction { vector, .. } => vector, + Load::Gravity { acceleration, .. } => acceleration, + Load::Pressure { magnitude, .. } => { + finite("pressure", *magnitude)?; + return Ok(()); + } + }; + for c in vector { + finite(self.kind(), *c)?; + } + if vector.iter().all(|c| *c == 0.0) { + return Err(Error::Invalid(format!("a {} vector must not be zero", self.kind()))); + } + Ok(()) + } +} + +impl Model { + pub fn loads(&self) -> &BTreeMap { + &self.loads + } + + pub fn load(&self, name: &str) -> Option<&Load> { + self.loads.get(name) + } + + fn check_load(&self, name: &str, load: &Load) -> Result<()> { + load.validate().map_err(|e| context("load", name, e))?; + if let Some(s) = load.set() { + require(&self.sets, "set", s).map_err(|e| context("load", name, e))?; + } + if let Some(r) = load.reference_point() { + require(&self.reference_points, "reference point", r).map_err(|e| context("load", name, e))?; + } + Ok(()) + } + + pub fn add_load(&mut self, name: &str, load: Load) -> Result<()> { + self.check_load(name, &load)?; + insert_new(&mut self.loads, "load", name, load) + } + + pub fn replace_load(&mut self, name: &str, load: Load) -> Result<()> { + self.check_load(name, &load)?; + replace_existing(&mut self.loads, "load", name, load) + } + + /// Remove a load that no step lists. + pub fn remove_load(&mut self, name: &str) -> Result<()> { + let users = self.steps.iter().filter(|s| s.loads.contains(name)).map(|s| format!("step `{}`", s.name)).collect(); + not_in_use("load", name, users)?; + remove_existing(&mut self.loads, "load", name).map(|_| ()) + } +} diff --git a/crates/model/src/materials.rs b/crates/model/src/materials.rs new file mode 100644 index 0000000..a4e850a --- /dev/null +++ b/crates/model/src/materials.rs @@ -0,0 +1,153 @@ +//! Materials: linear elasticity and mass density. + +use serde::{Deserialize, Serialize}; + +use crate::{Error, Model, Result, context, insert_new, not_in_use, positive, remove_existing, replace_existing}; + +/// Names of the engineering constants' shear moduli, in order. +pub(crate) const SHEAR_MODULI: [&str; 3] = ["G12", "G13", "G23"]; + +/// Linear elastic behaviour. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] +pub enum Elastic { + /// Young's modulus and Poisson's ratio. + Isotropic { young: f64, poisson: f64 }, + /// Orthotropic, by engineering constants in the material axes 1, 2, 3 (the global X, Y, Z + /// until material orientations land): Young's moduli `[E1, E2, E3]`, Poisson's ratios + /// `[ν12, ν13, ν23]` (ν_ij is the contraction in j under stress in i) and shear moduli + /// `[G12, G13, G23]`. + EngineeringConstants { young: [f64; 3], poisson: [f64; 3], shear: [f64; 3] }, +} + +impl Elastic { + /// Check the constants describe a stable material (positive-definite stiffness). + pub fn validate(&self) -> Result<()> { + match *self { + Elastic::Isotropic { young, poisson } => { + positive("Young's modulus", young)?; + if !(poisson > -1.0 && poisson < 0.5) { + return Err(Error::Invalid(format!("Poisson's ratio must be above -1 and below 0.5, got {poisson}"))); + } + } + Elastic::EngineeringConstants { young, poisson, shear } => { + for (i, e) in young.iter().enumerate() { + positive(&format!("E{}", i + 1), *e)?; + } + for (what, g) in SHEAR_MODULI.iter().zip(shear) { + positive(what, g)?; + } + let [e1, e2, e3] = young; + let [n12, n13, n23] = poisson; + if poisson.iter().any(|n| !n.is_finite()) { + return Err(Error::Invalid("Poisson's ratios must be finite numbers".into())); + } + // Reciprocal ratios ν_ji = ν_ij E_j / E_i; stability needs |ν_ij| < sqrt(E_i/E_j) + // and a positive determinant of the compliance's normal block. + let (n21, n31, n32) = (n12 * e2 / e1, n13 * e3 / e1, n23 * e3 / e2); + for (name, nij, ei, ej) in [("ν12", n12, e1, e2), ("ν13", n13, e1, e3), ("ν23", n23, e2, e3)] { + if !(nij.abs() < (ei / ej).sqrt()) { + return Err(Error::Invalid(format!("{name} = {nij} is too large for the moduli (must be below {})", (ei / ej).sqrt()))); + } + } + let det = 1.0 - n12 * n21 - n23 * n32 - n31 * n13 - 2.0 * n21 * n32 * n13; + if !(det > 0.0) { + return Err(Error::Invalid("the Poisson's ratios make the material unstable (stiffness not positive definite)".into())); + } + } + } + Ok(()) + } +} + +/// A material. +#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Material { + #[serde(default, skip_serializing_if = "Option::is_none")] + pub elastic: Option, + /// Mass per volume, needed by frequency steps and gravity. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub density: Option, +} + +impl Material { + pub fn validate(&self) -> Result<()> { + if self.elastic.is_none() && self.density.is_none() { + return Err(Error::Invalid("a material needs elastic constants or a density".into())); + } + if let Some(e) = &self.elastic { + e.validate()?; + } + if let Some(d) = self.density { + positive("density", d)?; + } + Ok(()) + } +} + +impl Model { + pub fn materials(&self) -> &std::collections::BTreeMap { + &self.materials + } + + pub fn material(&self, name: &str) -> Option<&Material> { + self.materials.get(name) + } + + pub fn add_material(&mut self, name: &str, material: Material) -> Result<()> { + material.validate().map_err(|e| context("material", name, e))?; + insert_new(&mut self.materials, "material", name, material) + } + + pub fn replace_material(&mut self, name: &str, material: Material) -> Result<()> { + material.validate().map_err(|e| context("material", name, e))?; + replace_existing(&mut self.materials, "material", name, material) + } + + pub fn remove_material(&mut self, name: &str) -> Result<()> { + let users = self.sections.iter().filter(|(_, s)| s.material == name).map(|(n, _)| format!("section `{n}`")).collect(); + not_in_use("material", name, users)?; + remove_existing(&mut self.materials, "material", name).map(|_| ()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn iso(young: f64, poisson: f64) -> Material { + Material { elastic: Some(Elastic::Isotropic { young, poisson }), density: None } + } + + #[test] + fn materials_validate_and_name_the_bad_constant() { + let mut m = Model::default(); + m.add_material("Steel", Material { density: Some(7.85e-9), ..iso(2.1e5, 0.3) }).unwrap(); + assert!(m.add_material("Steel", iso(1.0, 0.3)).is_err(), "duplicate"); + for bad in [iso(0.0, 0.3), iso(-1.0, 0.3), iso(f64::NAN, 0.3), iso(1.0, 0.5), iso(1.0, -1.0), iso(1.0, f64::INFINITY), Material::default()] { + let e = m.add_material("Bad", bad).unwrap_err(); + assert!(e.to_string().contains("material `Bad`"), "{e}"); + } + assert!(m.add_material("Light", Material { elastic: None, density: Some(-1.0) }).is_err()); + assert!(m.add_material("1st", iso(1.0, 0.3)).is_err(), "names start with a letter"); + assert!(m.add_material("has space", iso(1.0, 0.3)).is_err()); + m.replace_material("Steel", iso(2.0e5, 0.29)).unwrap(); + assert!(m.replace_material("Nope", iso(2.0e5, 0.29)).is_err()); + assert_eq!(m.materials().len(), 1); + } + + #[test] + fn orthotropic_stability() { + let ortho = |poisson| Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson, shear: [5e3, 5e3, 3.8e3] }; + ortho([0.3, 0.3, 0.4]).validate().unwrap(); + // ν21 = ν12 E2/E1 is small; ν12 itself may exceed 0.5 in a stiff direction. + ortho([0.9, 0.3, 0.4]).validate().unwrap(); + assert!(ortho([4.0, 0.3, 0.4]).validate().is_err(), "|ν12| < sqrt(E1/E2) = 3.74"); + assert!(ortho([0.3, 0.3, 1.0]).validate().is_err()); + assert!(ortho([3.5, 3.5, 0.9]).validate().is_err(), "determinant"); + assert!(ortho([0.3, f64::NAN, 0.4]).validate().is_err()); + let zero_shear = Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [1.0, 0.0, 1.0] }; + assert!(zero_shear.validate().unwrap_err().to_string().contains("G13")); + } +} diff --git a/crates/model/src/mesh.rs b/crates/model/src/mesh.rs new file mode 100644 index 0000000..f7a3035 --- /dev/null +++ b/crates/model/src/mesh.rs @@ -0,0 +1,400 @@ +//! Nodes and elements. + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::{Error, MAX_COORD, MAX_ELEMENTS, MAX_LABEL, MAX_NODES, Model, Result}; + +/// A node: a point in model space. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +pub struct Node { + pub xyz: [f64; 3], +} + +/// Element topology. Node order is corners first, then mid-side nodes, with positive volume for +/// a right-handed ordering: +/// - tet10: mid-side nodes 5–10 on edges 1-2, 2-3, 3-1, 1-4, 2-4, 3-4; +/// - wedge6: triangle 1-2-3, then 4-5-6 on the side its right-hand normal points to; +/// - hex8/hex20: face 1-2-3-4, then 5-6-7-8 above it; hex20's mid-side nodes 9–12 lie on the +/// edges of 1-2-3-4, 13–16 on those of 5-6-7-8, and 17–20 on edges 1-5, 2-6, 3-7, 4-8. +/// +/// This is the Abaqus (and VTK) order. Formats that order nodes differently, such as Nastran's +/// CHEXA with 20 nodes (vertical edges before the top face's edges), are reordered by their deck +/// readers and writers; the model never stores another order. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum ElementKind { + /// Two-node axial member. + Rod2, + /// Two-node beam. + Beam2, + /// Three-node triangle (shell or membrane). + Tri3, + /// Six-node triangle. + Tri6, + /// Four-node quadrilateral. + Quad4, + /// Eight-node quadrilateral. + Quad8, + /// Four-node tetrahedron. + Tet4, + /// Ten-node tetrahedron. + Tet10, + /// Six-node wedge (pentahedron). + Wedge6, + /// Eight-node hexahedron. + Hex8, + /// Twenty-node hexahedron. + Hex20, +} + +impl ElementKind { + pub const ALL: [ElementKind; 11] = [ + ElementKind::Rod2, + ElementKind::Beam2, + ElementKind::Tri3, + ElementKind::Tri6, + ElementKind::Quad4, + ElementKind::Quad8, + ElementKind::Tet4, + ElementKind::Tet10, + ElementKind::Wedge6, + ElementKind::Hex8, + ElementKind::Hex20, + ]; + + pub fn node_count(self) -> usize { + match self { + ElementKind::Rod2 | ElementKind::Beam2 => 2, + ElementKind::Tri3 => 3, + ElementKind::Quad4 | ElementKind::Tet4 => 4, + ElementKind::Tri6 | ElementKind::Wedge6 => 6, + ElementKind::Quad8 | ElementKind::Hex8 => 8, + ElementKind::Tet10 => 10, + ElementKind::Hex20 => 20, + } + } + + /// Topological dimension: 1 for line elements, 2 for shells, 3 for solids. + pub fn dimension(self) -> u8 { + match self { + ElementKind::Rod2 | ElementKind::Beam2 => 1, + ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => 2, + ElementKind::Tet4 | ElementKind::Tet10 | ElementKind::Wedge6 | ElementKind::Hex8 | ElementKind::Hex20 => 3, + } + } + + /// The name used in commands and `.modelcraft` files (`hex8`, `tet10`, …). + pub fn name(self) -> &'static str { + match self { + ElementKind::Rod2 => "rod2", + ElementKind::Beam2 => "beam2", + ElementKind::Tri3 => "tri3", + ElementKind::Tri6 => "tri6", + ElementKind::Quad4 => "quad4", + ElementKind::Quad8 => "quad8", + ElementKind::Tet4 => "tet4", + ElementKind::Tet10 => "tet10", + ElementKind::Wedge6 => "wedge6", + ElementKind::Hex8 => "hex8", + ElementKind::Hex20 => "hex20", + } + } + + pub fn from_name(name: &str) -> Option { + ElementKind::ALL.into_iter().find(|k| k.name() == name) + } + + /// The faces, numbered from 1 in this order, as positions in the element's node list: the + /// corners, then (quadratic kinds) the mid-side nodes in the same cyclic order. + /// + /// - Solids: each face's corners run counter-clockwise seen from outside, so the right-hand + /// normal points out of the element. Hexahedra: 1 is the face of nodes 1-2-3-4, 2 the face + /// of 5-6-7-8, then the four sides along edges 1-2, 2-3, 3-4 and 4-1. Tetrahedra: 1 is + /// 1-2-3, then the faces opposite node 3, node 1 and node 2. Wedges: the two triangles, then + /// the sides along edges 1-2, 2-3 and 3-1. + /// - Shells: face 1 is the side the element's right-hand normal points to, face 2 the other. + /// - Line elements have none. + pub fn faces(self) -> &'static [&'static [usize]] { + match self { + ElementKind::Rod2 | ElementKind::Beam2 => &[], + ElementKind::Tri3 => &[&[0, 1, 2], &[0, 2, 1]], + ElementKind::Tri6 => &[&[0, 1, 2, 3, 4, 5], &[0, 2, 1, 5, 4, 3]], + ElementKind::Quad4 => &[&[0, 1, 2, 3], &[0, 3, 2, 1]], + ElementKind::Quad8 => &[&[0, 1, 2, 3, 4, 5, 6, 7], &[0, 3, 2, 1, 7, 6, 5, 4]], + ElementKind::Tet4 => &[&[0, 2, 1], &[0, 1, 3], &[1, 2, 3], &[0, 3, 2]], + ElementKind::Tet10 => &[&[0, 2, 1, 6, 5, 4], &[0, 1, 3, 4, 8, 7], &[1, 2, 3, 5, 9, 8], &[0, 3, 2, 7, 9, 6]], + ElementKind::Wedge6 => &[&[0, 2, 1], &[3, 4, 5], &[0, 1, 4, 3], &[1, 2, 5, 4], &[2, 0, 3, 5]], + ElementKind::Hex8 => &[&[0, 3, 2, 1], &[4, 5, 6, 7], &[0, 1, 5, 4], &[1, 2, 6, 5], &[2, 3, 7, 6], &[3, 0, 4, 7]], + ElementKind::Hex20 => &[ + &[0, 3, 2, 1, 11, 10, 9, 8], + &[4, 5, 6, 7, 12, 13, 14, 15], + &[0, 1, 5, 4, 8, 17, 12, 16], + &[1, 2, 6, 5, 9, 18, 13, 17], + &[2, 3, 7, 6, 10, 19, 14, 18], + &[3, 0, 4, 7, 11, 16, 15, 19], + ], + } + } + + /// The edges, numbered from 1 in this order, as positions in the element's node list: the two + /// corners, then the mid-side node for quadratic kinds. The order follows the mid-side nodes + /// (see [`ElementKind`]); linear kinds use the same order. + pub fn edges(self) -> &'static [&'static [usize]] { + match self { + ElementKind::Rod2 | ElementKind::Beam2 => &[&[0, 1]], + ElementKind::Tri3 => &[&[0, 1], &[1, 2], &[2, 0]], + ElementKind::Tri6 => &[&[0, 1, 3], &[1, 2, 4], &[2, 0, 5]], + ElementKind::Quad4 => &[&[0, 1], &[1, 2], &[2, 3], &[3, 0]], + ElementKind::Quad8 => &[&[0, 1, 4], &[1, 2, 5], &[2, 3, 6], &[3, 0, 7]], + ElementKind::Tet4 => &[&[0, 1], &[1, 2], &[2, 0], &[0, 3], &[1, 3], &[2, 3]], + ElementKind::Tet10 => &[&[0, 1, 4], &[1, 2, 5], &[2, 0, 6], &[0, 3, 7], &[1, 3, 8], &[2, 3, 9]], + ElementKind::Wedge6 => &[&[0, 1], &[1, 2], &[2, 0], &[3, 4], &[4, 5], &[5, 3], &[0, 3], &[1, 4], &[2, 5]], + ElementKind::Hex8 => &[&[0, 1], &[1, 2], &[2, 3], &[3, 0], &[4, 5], &[5, 6], &[6, 7], &[7, 4], &[0, 4], &[1, 5], &[2, 6], &[3, 7]], + ElementKind::Hex20 => &[ + &[0, 1, 8], + &[1, 2, 9], + &[2, 3, 10], + &[3, 0, 11], + &[4, 5, 12], + &[5, 6, 13], + &[6, 7, 14], + &[7, 4, 15], + &[0, 4, 16], + &[1, 5, 17], + &[2, 6, 18], + &[3, 7, 19], + ], + } + } +} + +/// One face of an element, a member of a surface: the element's label and its face number +/// (1-based, [`ElementKind::faces`]). +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct ElementFace { + pub element: u32, + pub face: u8, +} + +/// An element: a topology and its node labels. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +pub struct Element { + pub kind: ElementKind, + pub nodes: Vec, +} + +/// Check that every coordinate is a finite number within ±[`MAX_COORD`] (the rule for nodes, +/// shared by everything that takes coordinates from outside). +pub fn check_coordinates(xyz: &[f64]) -> Result<()> { + match xyz.iter().find(|c| !(c.abs() <= MAX_COORD)) { + Some(c) => Err(Error::Invalid(format!("coordinate {c} is not a finite number within ±{MAX_COORD}"))), + None => Ok(()), + } +} + +pub(crate) fn check_label(what: &str, label: u32) -> Result<()> { + if label == 0 || label > MAX_LABEL { + return Err(Error::Invalid(format!("{what} label {label} is outside 1..={MAX_LABEL}"))); + } + Ok(()) +} + +impl Model { + pub fn nodes(&self) -> &BTreeMap { + &self.nodes + } + + pub fn elements(&self) -> &BTreeMap { + &self.elements + } + + pub fn node(&self, label: u32) -> Option<&Node> { + self.nodes.get(&label) + } + + pub fn element(&self, label: u32) -> Option<&Element> { + self.elements.get(&label) + } + + /// Add a node, or move an existing one (decks redefine nodes the same way). + pub fn set_node(&mut self, label: u32, xyz: [f64; 3]) -> Result<()> { + check_label("node", label)?; + check_coordinates(&xyz).map_err(|e| match e { + Error::Invalid(m) => Error::Invalid(format!("node {label}: {m}")), + e => e, + })?; + if !self.nodes.contains_key(&label) && self.nodes.len() >= MAX_NODES { + return Err(Error::Limit(format!("a model holds at most {MAX_NODES} nodes"))); + } + self.nodes.insert(label, Node { xyz }); + Ok(()) + } + + /// Add an element, or replace an existing one. Its nodes must exist and be distinct. + pub fn set_element(&mut self, label: u32, kind: ElementKind, nodes: Vec) -> Result<()> { + check_label("element", label)?; + if nodes.len() != kind.node_count() { + return Err(Error::Invalid(format!("element {label}: {} needs {} nodes, got {}", kind.name(), kind.node_count(), nodes.len()))); + } + for (i, n) in nodes.iter().enumerate() { + if !self.nodes.contains_key(n) { + return Err(Error::NotFound(format!("element {label}: node {n}"))); + } + if nodes.get(..i).is_some_and(|before| before.contains(n)) { + return Err(Error::Invalid(format!("element {label}: node {n} appears twice"))); + } + } + if !self.elements.contains_key(&label) && self.elements.len() >= MAX_ELEMENTS { + return Err(Error::Limit(format!("a model holds at most {MAX_ELEMENTS} elements"))); + } + // A replaced element's face and edge numbers no longer mean what its tags recorded. + self.classification.forget_element(label); + self.elements.insert(label, Element { kind, nodes }); + Ok(()) + } + + /// Remove a node that no element uses. + pub fn remove_node(&mut self, label: u32) -> Result<()> { + if let Some((e, _)) = self.elements.iter().find(|(_, e)| e.nodes.contains(&label)) { + return Err(Error::Invalid(format!("node {label} is used by element {e}"))); + } + self.nodes.remove(&label).ok_or_else(|| Error::NotFound(format!("node {label}")))?; + self.classification.forget_node(label); + Ok(()) + } + + pub fn remove_element(&mut self, label: u32) -> Result<()> { + self.elements.remove(&label).ok_or_else(|| Error::NotFound(format!("element {label}")))?; + self.classification.forget_element(label); + Ok(()) + } + + /// Element counts by kind, in [`ElementKind::ALL`] order, kinds with none left out. + pub fn element_counts(&self) -> Vec<(ElementKind, usize)> { + ElementKind::ALL.into_iter().map(|k| (k, self.elements.values().filter(|e| e.kind == k).count())).filter(|(_, n)| *n > 0).collect() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn cube() -> Model { + let mut m = Model::default(); + let corners = [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]]; + for (i, c) in corners.into_iter().enumerate() { + m.set_node(i as u32 + 1, c).unwrap(); + } + m.set_element(1, ElementKind::Hex8, (1..=8).collect()).unwrap(); + m + } + + #[test] + fn builds_a_hex() { + let m = cube(); + assert_eq!(m.nodes().len(), 8); + assert_eq!(m.element_counts(), vec![(ElementKind::Hex8, 1)]); + } + + #[test] + fn rejects_bad_labels_and_coordinates() { + let mut m = Model::default(); + assert!(m.set_node(0, [0.; 3]).is_err()); + assert!(m.set_node(MAX_LABEL + 1, [0.; 3]).is_err()); + assert!(m.set_node(1, [f64::NAN, 0., 0.]).is_err()); + assert!(m.set_node(1, [f64::INFINITY, 0., 0.]).is_err()); + assert!(m.set_node(1, [2e9, 0., 0.]).is_err()); + assert!(m.nodes().is_empty()); + } + + #[test] + fn rejects_bad_connectivity() { + let mut m = cube(); + assert!(matches!(m.set_element(2, ElementKind::Tet4, vec![1, 2, 3]), Err(Error::Invalid(_)))); + assert!(matches!(m.set_element(2, ElementKind::Tet4, vec![1, 2, 3, 99]), Err(Error::NotFound(_)))); + assert!(matches!(m.set_element(2, ElementKind::Tet4, vec![1, 2, 3, 3]), Err(Error::Invalid(_)))); + assert_eq!(m.elements().len(), 1); + } + + #[test] + fn used_nodes_cannot_be_removed() { + let mut m = cube(); + assert!(m.remove_node(1).is_err()); + m.remove_element(1).unwrap(); + m.remove_node(1).unwrap(); + assert!(m.remove_node(1).is_err()); + } + + /// Reference coordinates with positive volume for each kind (mid-side nodes at edge midpoints). + fn reference(kind: ElementKind) -> Vec<[f64; 3]> { + let corners: Vec<[f64; 3]> = match kind { + ElementKind::Rod2 | ElementKind::Beam2 => vec![[0., 0., 0.], [1., 0., 0.]], + ElementKind::Tri3 | ElementKind::Tri6 => 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.]], + ElementKind::Hex8 | ElementKind::Hex20 => { + 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(), [f64::NAN; 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); + } + } + xyz + } + + #[test] + fn edges_cover_every_node_and_faces_point_out() { + for kind in ElementKind::ALL { + let xyz = reference(kind); + assert!(xyz.iter().flatten().all(|c| c.is_finite()), "{kind:?}: every mid-side node is on an edge"); + let centroid = [0, 1, 2].map(|i| xyz.iter().map(|p| p[i]).sum::() / xyz.len() as f64); + for (f, face) in kind.faces().iter().enumerate() { + let n = if face.len() > 4 { face.len() / 2 } else { face.len() }; + let c: Vec<[f64; 3]> = face[..n].iter().map(|&i| xyz[i]).collect(); + // Newell normal of the corner polygon, against the direction out of the centroid. + let mut normal = [0.0; 3]; + for k in 0..n { + let (p, q) = (c[k], c[(k + 1) % n]); + normal = [normal[0] + p[1] * q[2] - p[2] * q[1], normal[1] + p[2] * q[0] - p[0] * q[2], normal[2] + p[0] * q[1] - p[1] * q[0]]; + } + let out = [0, 1, 2].map(|i| c[0][i] - centroid[i]); + if kind.dimension() == 3 { + assert!(normal[0] * out[0] + normal[1] * out[1] + normal[2] * out[2] > 0.0, "{kind:?} face {} points in", f + 1); + } + // A quadratic face's mid-side nodes sit on its own corner edges, in cyclic order. + for (k, &m) in face[n..].iter().enumerate() { + let (a, b) = (face[k], face[(k + 1) % n]); + assert_eq!(xyz[m], [0, 1, 2].map(|i| (xyz[a][i] + xyz[b][i]) / 2.0), "{kind:?} face {} mid-side {}", f + 1, k + 1); + } + } + if kind.dimension() == 2 { + let [p, q] = kind.faces() else { panic!("{kind:?}: two sides") }; + assert_eq!(p.len(), kind.node_count()); + let mut sorted = q.to_vec(); + sorted.sort_unstable(); + assert_eq!(sorted, (0..kind.node_count()).collect::>()); + } + } + assert_eq!(ElementKind::Hex8.faces().len(), 6); + assert_eq!(ElementKind::Tet4.faces().len(), 4); + assert_eq!(ElementKind::Wedge6.faces().len(), 5); + assert_eq!(ElementKind::Hex20.edges().len(), 12); + } + + #[test] + fn kind_names_round_trip() { + for k in ElementKind::ALL { + assert_eq!(ElementKind::from_name(k.name()), Some(k)); + assert_eq!(serde_json::to_value(k).unwrap(), serde_json::json!(k.name())); + } + assert_eq!(ElementKind::from_name("C3D8"), None); + } +} diff --git a/crates/model/src/nodal.rs b/crates/model/src/nodal.rs new file mode 100644 index 0000000..613c6e8 --- /dev/null +++ b/crates/model/src/nodal.rs @@ -0,0 +1,183 @@ +//! Consistent nodal loads: turning a force, moment or traction on a set into the nodal values a +//! deck or the solver applies, the same way for every consumer. +//! +//! A load over a surface is spread by the element shape functions: node `i` of a face takes +//! `∫ N_i dA` of a unit traction, integrated by Gauss quadrature on the face's own (possibly +//! curved, quadratic) geometry. The weights of a face sum to its area, so a total force split by +//! them keeps its total on any mesh of the same surface. Quadratic faces give the textbook +//! values (corner 0 and mid-side A/3 on a flat six-node triangle; corner −A/12 and mid-side A/3 on +//! a flat eight-node quadrilateral). + +use std::collections::BTreeMap; + +use crate::{ElementFace, Error, Model, Result}; + +/// Shape functions and their derivatives at one point of a reference face. +struct Shape { + n: Vec, + dxi: Vec, + deta: Vec, +} + +/// Triangle: barycentric L0 = 1 − ξ − η, L1 = ξ, L2 = η. +fn triangle(nodes: usize, xi: f64, eta: f64) -> Shape { + let l = [1.0 - xi - eta, xi, eta]; + let dl = [(-1.0, -1.0), (1.0, 0.0), (0.0, 1.0)]; + let mut s = Shape { n: Vec::with_capacity(nodes), dxi: Vec::with_capacity(nodes), deta: Vec::with_capacity(nodes) }; + if nodes == 3 { + for (li, (dx, de)) in l.iter().zip(dl) { + s.n.push(*li); + s.dxi.push(dx); + s.deta.push(de); + } + return s; + } + for (li, (dx, de)) in l.iter().zip(dl) { + s.n.push(li * (2.0 * li - 1.0)); + s.dxi.push((4.0 * li - 1.0) * dx); + s.deta.push((4.0 * li - 1.0) * de); + } + for (la, lb, (dxa, dea), (dxb, deb)) in [(l[0], l[1], dl[0], dl[1]), (l[1], l[2], dl[1], dl[2]), (l[2], l[0], dl[2], dl[0])] { + s.n.push(4.0 * la * lb); + s.dxi.push(4.0 * (dxa * lb + la * dxb)); + s.deta.push(4.0 * (dea * lb + la * deb)); + } + s +} + +/// Quadrilateral on [−1, 1]², corners (−1,−1), (1,−1), (1,1), (−1,1). +fn quadrilateral(nodes: usize, xi: f64, eta: f64) -> Shape { + const CORNERS: [(f64, f64); 4] = [(-1.0, -1.0), (1.0, -1.0), (1.0, 1.0), (-1.0, 1.0)]; + let mut s = Shape { n: Vec::with_capacity(nodes), dxi: Vec::with_capacity(nodes), deta: Vec::with_capacity(nodes) }; + for (xi_i, eta_i) in CORNERS { + let (a, b) = (1.0 + xi * xi_i, 1.0 + eta * eta_i); + if nodes == 4 { + s.n.push(a * b / 4.0); + s.dxi.push(xi_i * b / 4.0); + s.deta.push(eta_i * a / 4.0); + } else { + let c = xi * xi_i + eta * eta_i - 1.0; + s.n.push(a * b * c / 4.0); + s.dxi.push(xi_i * b * (2.0 * xi * xi_i + eta * eta_i) / 4.0); + s.deta.push(eta_i * a * (xi * xi_i + 2.0 * eta * eta_i) / 4.0); + } + } + if nodes == 8 { + // Mid-sides of edges 1-2 (η = −1), 2-3 (ξ = 1), 3-4 (η = 1), 4-1 (ξ = −1). + for (xi_i, eta_i) in [(0.0, -1.0), (1.0, 0.0), (0.0, 1.0), (-1.0, 0.0)] { + if xi_i == 0.0 { + s.n.push((1.0 - xi * xi) * (1.0 + eta * eta_i) / 2.0); + s.dxi.push(-xi * (1.0 + eta * eta_i)); + s.deta.push(eta_i * (1.0 - xi * xi) / 2.0); + } else { + s.n.push((1.0 + xi * xi_i) * (1.0 - eta * eta) / 2.0); + s.dxi.push(xi_i * (1.0 - eta * eta) / 2.0); + s.deta.push(-eta * (1.0 + xi * xi_i)); + } + } + } + s +} + +/// Gauss points `(ξ, η, weight)`: a degree-4 six-point rule on the triangle (Dunavant 1985) and +/// 3 × 3 Gauss–Legendre on the quadrilateral; both integrate the quadratic faces' weights exactly +/// when the face is flat. +fn rule(corners: usize) -> Vec<(f64, f64, f64)> { + if corners == 3 { + let (a1, w1) = (0.445_948_490_915_965, 0.223_381_589_678_011 / 2.0); + let (a2, w2) = (0.091_576_213_509_771, 0.109_951_743_655_322 / 2.0); + let (b1, b2) = (1.0 - 2.0 * a1, 1.0 - 2.0 * a2); + vec![(a1, a1, w1), (a1, b1, w1), (b1, a1, w1), (a2, a2, w2), (a2, b2, w2), (b2, a2, w2)] + } else { + let g = (0.6f64).sqrt(); + let pts = [(-g, 5.0 / 9.0), (0.0, 8.0 / 9.0), (g, 5.0 / 9.0)]; + pts.iter().flat_map(|&(x, wx)| pts.iter().map(move |&(y, wy)| (x, y, wx * wy))).collect() + } +} + +impl Model { + /// The nodal weights `∫ N_i dA` of a set of element faces, summed per node. Their sum is the + /// faces' area. + pub fn surface_weights(&self, faces: &[ElementFace]) -> Result> { + let mut out: BTreeMap = BTreeMap::new(); + for &f in faces { + let labels = self.face_nodes(f).ok_or_else(|| Error::NotFound(format!("element {} face {}", f.element, f.face)))?; + let xyz: Vec<[f64; 3]> = + labels.iter().map(|n| self.node(*n).map(|p| p.xyz).ok_or_else(|| Error::NotFound(format!("node {n}")))).collect::>()?; + let (corners, shape): (usize, fn(usize, f64, f64) -> Shape) = match labels.len() { + 3 | 6 => (3, triangle), + 4 | 8 => (4, quadrilateral), + n => return Err(Error::Invalid(format!("element {} face {} has {n} nodes, not a surface face", f.element, f.face))), + }; + for (xi, eta, w) in rule(corners) { + let s = shape(labels.len(), xi, eta); + let mut tx = [0.0; 3]; + let mut te = [0.0; 3]; + for ((p, dx), de) in xyz.iter().zip(&s.dxi).zip(&s.deta) { + for k in 0..3 { + tx[k] += dx * p[k]; + te[k] += de * p[k]; + } + } + let c = [tx[1] * te[2] - tx[2] * te[1], tx[2] * te[0] - tx[0] * te[2], tx[0] * te[1] - tx[1] * te[0]]; + let da = (c[0] * c[0] + c[1] * c[1] + c[2] * c[2]).sqrt() * w; + for (label, n) in labels.iter().zip(&s.n) { + *out.entry(*label).or_default() += n * da; + } + } + } + Ok(out) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::ElementKind; + + /// A flat 2 × 3 face of a hex8 or hex20 (face 1, z = 0), and its weights. + fn face_weights(kind: ElementKind) -> Vec { + let mut m = Model::default(); + let corners = [[0., 0., 0.], [2., 0., 0.], [2., 3., 0.], [0., 3., 0.], [0., 0., 1.], [2., 0., 1.], [2., 3., 1.], [0., 3., 1.]]; + let mut xyz = corners.to_vec(); + for e in kind.edges() { + if let [a, b, _] = **e { + xyz.push([0, 1, 2].map(|i| (corners[a][i] + corners[b][i]) / 2.0)); + } + } + for (i, p) in xyz.iter().enumerate() { + m.set_node(i as u32 + 1, *p).unwrap(); + } + m.set_element(1, kind, (1..=kind.node_count() as u32).collect()).unwrap(); + let w = m.surface_weights(&[ElementFace { element: 1, face: 1 }]).unwrap(); + m.face_nodes(ElementFace { element: 1, face: 1 }).unwrap().iter().map(|n| w[n]).collect() + } + + #[test] + fn weights_are_the_consistent_nodal_loads() { + let close = |a: &[f64], b: &[f64]| a.iter().zip(b).all(|(x, y)| (x - y).abs() < 1e-12); + assert!(close(&face_weights(ElementKind::Hex8), &[1.5; 4]), "{:?}", face_weights(ElementKind::Hex8)); + let quad8 = face_weights(ElementKind::Hex20); + assert!(close(&quad8, &[-0.5, -0.5, -0.5, -0.5, 2.0, 2.0, 2.0, 2.0]), "{quad8:?}"); + + // Triangles: a tet4 and a tet10 face (face 1, the right triangle of legs 1 in z = 0). + let mut m = Model::default(); + let p = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]; + let mut xyz = p.to_vec(); + for e in ElementKind::Tet10.edges() { + if let [a, b, _] = **e { + xyz.push([0, 1, 2].map(|i| (p[a][i] + p[b][i]) / 2.0)); + } + } + for (i, q) in xyz.iter().enumerate() { + m.set_node(i as u32 + 1, *q).unwrap(); + } + m.set_element(1, ElementKind::Tet4, vec![1, 2, 3, 4]).unwrap(); + m.set_element(2, ElementKind::Tet10, (1..=10).collect()).unwrap(); + let w = m.surface_weights(&[ElementFace { element: 1, face: 1 }]).unwrap(); + assert!(close(&w.values().copied().collect::>(), &[1.0 / 6.0; 3]), "{w:?}"); + let w = m.surface_weights(&[ElementFace { element: 2, face: 1 }]).unwrap(); + let tri6: Vec = m.face_nodes(ElementFace { element: 2, face: 1 }).unwrap().iter().map(|n| w[n]).collect(); + assert!(close(&tri6, &[0.0, 0.0, 0.0, 1.0 / 6.0, 1.0 / 6.0, 1.0 / 6.0]), "{tri6:?}"); + } +} diff --git a/crates/model/src/sections.rs b/crates/model/src/sections.rs new file mode 100644 index 0000000..3afc2f1 --- /dev/null +++ b/crates/model/src/sections.rs @@ -0,0 +1,306 @@ +//! Sections: what the elements of a set are made of, how thick they are, and how they are +//! integrated. +//! +//! A section is assigned to one set and carries the element formulation for its elements, so one +//! region has one formulation and every solver can express it: the writers map the element kind +//! plus the [`Formulation`] to each format's element names, property cards or formulation +//! numbers. + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::{Error, Model, Result, context, finite, insert_new, positive, remove_existing, replace_existing, require}; + +/// Most integration points through a shell's thickness. +pub const MAX_SHELL_POINTS: u32 = 15; + +/// A beam's cross-section shape. Dimensions are in the section's local axes: `width` along the +/// beam's orientation vector (local 1), `height` across it (local 2). +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] +pub enum BeamProfile { + Rectangle { + width: f64, + height: f64, + }, + Circle { + radius: f64, + }, + /// A tube: outer radius and wall thickness. + Pipe { + radius: f64, + thickness: f64, + }, + /// A doubly symmetric I: overall height, flange width, web and flange thicknesses. + I { + height: f64, + width: f64, + web: f64, + flange: f64, + }, + /// A rectangular tube with one wall thickness. + Box { + width: f64, + height: f64, + thickness: f64, + }, +} + +impl BeamProfile { + pub fn validate(&self) -> Result<()> { + match *self { + BeamProfile::Rectangle { width, height } => { + positive("width", width)?; + positive("height", height)?; + } + BeamProfile::Circle { radius } => { + positive("radius", radius)?; + } + BeamProfile::Pipe { radius, thickness } => { + positive("radius", radius)?; + positive("thickness", thickness)?; + if thickness > radius { + return Err(Error::Invalid(format!("a pipe's wall ({thickness}) cannot be thicker than its radius ({radius})"))); + } + } + BeamProfile::I { height, width, web, flange } => { + for (what, v) in [("height", height), ("width", width), ("web", web), ("flange", flange)] { + positive(what, v)?; + } + if web >= width || 2.0 * flange >= height { + return Err(Error::Invalid( + "an I profile's web must be narrower than its flanges, and its flanges thinner than half its height".into(), + )); + } + } + BeamProfile::Box { width, height, thickness } => { + for (what, v) in [("width", width), ("height", height), ("thickness", thickness)] { + positive(what, v)?; + } + if 2.0 * thickness >= width.min(height) { + return Err(Error::Invalid("a box's walls must be thinner than half its width and height".into())); + } + } + } + Ok(()) + } + + /// Cross-sectional area. + pub fn area(&self) -> f64 { + match *self { + BeamProfile::Rectangle { width, height } => width * height, + BeamProfile::Circle { radius } => std::f64::consts::PI * radius * radius, + BeamProfile::Pipe { radius, thickness } => std::f64::consts::PI * (radius * radius - (radius - thickness).powi(2)), + BeamProfile::I { height, width, web, flange } => 2.0 * width * flange + (height - 2.0 * flange) * web, + BeamProfile::Box { width, height, thickness } => width * height - (width - 2.0 * thickness) * (height - 2.0 * thickness), + } + } +} + +/// What a section describes. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] +pub enum SectionKind { + /// Continuum solid elements. + Solid {}, + /// Shell elements: thickness, the reference surface's offset from the mid-surface as a + /// fraction of the thickness (0 mid-surface, 0.5 on the side the normal points to, -0.5 on + /// the other), and integration points through the thickness. + Shell { + thickness: f64, + #[serde(default)] + offset: f64, + #[serde(default = "default_shell_points")] + points: u32, + }, + /// Beam elements with a library profile. `orientation` is a vector giving the profile's + /// local 1 direction (projected normal to each beam's axis); it must not be parallel to a + /// beam. + Beam { profile: BeamProfile, orientation: [f64; 3] }, + /// Axial members (rods, trusses) with a cross-sectional area. + Rod { area: f64 }, +} + +fn default_shell_points() -> u32 { + 5 +} + +impl SectionKind { + /// The element dimension this section applies to. + pub fn dimension(&self) -> u8 { + match self { + SectionKind::Solid {} => 3, + SectionKind::Shell { .. } => 2, + SectionKind::Beam { .. } | SectionKind::Rod { .. } => 1, + } + } + + /// `solid`, `shell`, `beam` or `rod`. + pub fn name(&self) -> &'static str { + match self { + SectionKind::Solid {} => "solid", + SectionKind::Shell { .. } => "shell", + SectionKind::Beam { .. } => "beam", + SectionKind::Rod { .. } => "rod", + } + } + + pub fn validate(&self) -> Result<()> { + match self { + SectionKind::Solid {} => {} + SectionKind::Shell { thickness, offset, points } => { + positive("thickness", *thickness)?; + if !(finite("offset", *offset)?.abs() <= 0.5) { + return Err(Error::Invalid(format!("a shell offset is a fraction of the thickness in -0.5..=0.5, got {offset}"))); + } + if !(1..=MAX_SHELL_POINTS).contains(points) { + return Err(Error::Invalid(format!( + "a shell has 1 to {MAX_SHELL_POINTS} integration points through its thickness, got {points}" + ))); + } + } + SectionKind::Beam { profile, orientation } => { + profile.validate()?; + for c in orientation { + finite("orientation", *c)?; + } + if orientation.iter().all(|c| *c == 0.0) { + return Err(Error::Invalid("a beam's orientation vector must not be zero".into())); + } + } + SectionKind::Rod { area } => { + positive("area", *area)?; + } + } + Ok(()) + } +} + +/// Number of integration points per element direction. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum Integration { + /// The full rule for the element's order. + #[default] + Full, + /// One order lower (cheaper, no shear locking, needs hourglass control). + Reduced, +} + +/// How the section's elements are formulated, independent of any solver's element names. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Formulation { + #[serde(default)] + pub integration: Integration, + /// Enhanced (incompatible) bending modes: linear hexahedra only. + #[serde(default)] + pub incompatible_modes: bool, + /// Mixed displacement–pressure formulation for nearly incompressible materials: solids only. + #[serde(default)] + pub hybrid: bool, +} + +impl Formulation { + /// Check the options against each other and against the section kind (the element kinds are + /// checked when the model is resolved). + pub fn validate_for(&self, kind: &SectionKind) -> Result<()> { + if self.incompatible_modes && self.integration == Integration::Reduced { + return Err(Error::Invalid("incompatible modes use full integration".into())); + } + if kind.dimension() != 3 && (self.incompatible_modes || self.hybrid) { + return Err(Error::Invalid(format!("incompatible modes and hybrid formulations are for solid sections, not {}", kind.name()))); + } + Ok(()) + } +} + +/// A section: a material and a kind, assigned to the elements of a set. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Section { + pub set: String, + pub material: String, + pub kind: SectionKind, + #[serde(default)] + pub formulation: Formulation, +} + +impl Section { + pub fn validate(&self) -> Result<()> { + self.kind.validate()?; + self.formulation.validate_for(&self.kind) + } +} + +impl Model { + pub fn sections(&self) -> &BTreeMap { + &self.sections + } + + pub fn section(&self, name: &str) -> Option<&Section> { + self.sections.get(name) + } + + fn check_section(&self, name: &str, section: &Section) -> Result<()> { + section.validate().map_err(|e| context("section", name, e))?; + require(&self.sets, "set", §ion.set).map_err(|e| context("section", name, e))?; + require(&self.materials, "material", §ion.material).map_err(|e| context("section", name, e)) + } + + pub fn add_section(&mut self, name: &str, section: Section) -> Result<()> { + self.check_section(name, §ion)?; + insert_new(&mut self.sections, "section", name, section) + } + + pub fn replace_section(&mut self, name: &str, section: Section) -> Result<()> { + self.check_section(name, §ion)?; + replace_existing(&mut self.sections, "section", name, section) + } + + pub fn remove_section(&mut self, name: &str) -> Result<()> { + remove_existing(&mut self.sections, "section", name).map(|_| ()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn profiles_and_formulations_validate() { + assert!(BeamProfile::Pipe { radius: 1.0, thickness: 1.5 }.validate().is_err()); + assert!(BeamProfile::I { height: 10.0, width: 5.0, web: 6.0, flange: 1.0 }.validate().is_err()); + assert!(BeamProfile::Box { width: 10.0, height: 4.0, thickness: 2.0 }.validate().is_err()); + assert!(BeamProfile::Rectangle { width: 1.0, height: f64::NAN }.validate().is_err()); + let i = BeamProfile::I { height: 10.0, width: 5.0, web: 1.0, flange: 1.0 }; + i.validate().unwrap(); + assert_eq!(i.area(), 2.0 * 5.0 + 8.0); + assert!((BeamProfile::Pipe { radius: 2.0, thickness: 2.0 }.area() - BeamProfile::Circle { radius: 2.0 }.area()).abs() < 1e-12); + + let shell = SectionKind::Shell { thickness: 2.0, offset: 0.0, points: 5 }; + assert!(Formulation { hybrid: true, ..Formulation::default() }.validate_for(&shell).is_err()); + assert!( + Formulation { incompatible_modes: true, integration: Integration::Reduced, hybrid: false }.validate_for(&SectionKind::Solid {}).is_err() + ); + Formulation { incompatible_modes: true, hybrid: true, ..Formulation::default() }.validate_for(&SectionKind::Solid {}).unwrap(); + for bad in [ + SectionKind::Shell { thickness: 2.0, offset: 0.6, points: 5 }, + SectionKind::Shell { thickness: 2.0, offset: 0.0, points: 0 }, + SectionKind::Shell { thickness: 0.0, offset: 0.0, points: 5 }, + SectionKind::Beam { profile: BeamProfile::Circle { radius: 1.0 }, orientation: [0.0; 3] }, + SectionKind::Rod { area: -1.0 }, + ] { + assert!(bad.validate().is_err(), "{bad:?}"); + } + let parsed: SectionKind = serde_json::from_value(serde_json::json!({"type": "shell", "thickness": 1.5})).unwrap(); + assert_eq!(parsed, SectionKind::Shell { thickness: 1.5, offset: 0.0, points: 5 }); + assert!( + serde_json::from_value::(serde_json::json!({"type": "shell", "thickness": 1.5, "thicknes": 2})).is_err(), + "typos are errors" + ); + assert!(serde_json::from_value::(serde_json::json!({"type": "solid", "x": 1})).is_err()); + assert!(serde_json::from_value::(serde_json::json!({"type": "membrane"})).is_err()); + } +} diff --git a/crates/model/src/sets.rs b/crates/model/src/sets.rs new file mode 100644 index 0000000..af62d71 --- /dev/null +++ b/crates/model/src/sets.rs @@ -0,0 +1,253 @@ +//! Named sets: regions of the model that sections, couplings, boundary conditions and loads +//! apply to. +//! +//! A set is defined by geometry names, explicit labels or both; it is resolved against the +//! current mesh into nodes, elements and element faces ([`ResolvedSet`]). Defining sets by +//! geometry is the norm: the definition outlives the mesh, so a remesh (or an edit of the design +//! followed by a remesh) keeps every set, section and load in place. Explicit labels are for +//! meshes with no geometry, such as one read from a deck. +//! +//! What each member contributes: +//! - a body: the elements filling it, their nodes, and the element faces on the body's faces; +//! - a face: the element faces on it, the shell elements filling it, and their nodes; +//! - an edge: the nodes of element edges on it, the line elements filling it, and their nodes; +//! - a vertex: the node on it; +//! - node labels: those nodes; element labels: those elements and their nodes; element faces: +//! those faces and their nodes. + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::check::ResolvedSet; +use crate::mesh::check_label; +use crate::{ + ElementFace, Error, GeometryRef, MAX_ELEMENTS, MAX_NODES, Model, Result, context, insert_new, not_in_use, remove_existing, replace_existing, +}; + +/// Most geometry references in one set. +pub const MAX_SET_REFERENCES: usize = 100_000; + +/// A set's definition. Members are kept sorted and without repeats. +#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Set { + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub geometry: Vec, + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub nodes: Vec, + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub elements: Vec, + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub faces: Vec, +} + +impl Set { + /// A set of geometry entities. + pub fn of_geometry(geometry: Vec) -> Set { + Set { geometry, ..Set::default() } + } + + pub fn is_empty(&self) -> bool { + self.geometry.is_empty() && self.nodes.is_empty() && self.elements.is_empty() && self.faces.is_empty() + } + + /// Check the definition's size and labels, and sort its members. + pub fn normalized(mut self) -> Result { + if self.is_empty() { + return Err(Error::Invalid("a set needs geometry, nodes, elements or element faces".into())); + } + if self.geometry.len() > MAX_SET_REFERENCES { + return Err(Error::Limit(format!("a set names at most {MAX_SET_REFERENCES} geometry entities"))); + } + if self.nodes.len() > MAX_NODES || self.elements.len() > MAX_ELEMENTS || self.faces.len() > MAX_ELEMENTS { + return Err(Error::Limit("a set lists at most as many members as a model holds".into())); + } + for &l in &self.nodes { + check_label("node", l)?; + } + for &l in self.elements.iter().chain(self.faces.iter().map(|f| &f.element)) { + check_label("element", l)?; + } + if let Some(f) = self.faces.iter().find(|f| f.face == 0) { + return Err(Error::Invalid(format!("element {} face 0: faces are numbered from 1", f.element))); + } + for v in [&mut self.nodes, &mut self.elements] { + v.sort_unstable(); + v.dedup(); + } + self.faces.sort_unstable(); + self.faces.dedup(); + self.geometry.sort(); + self.geometry.dedup(); + Ok(self) + } +} + +impl Model { + pub fn sets(&self) -> &BTreeMap { + &self.sets + } + + pub fn set(&self, name: &str) -> Option<&Set> { + self.sets.get(name) + } + + pub fn add_set(&mut self, name: &str, set: Set) -> Result<()> { + let set = set.normalized().map_err(|e| context("set", name, e))?; + insert_new(&mut self.sets, "set", name, set) + } + + pub fn replace_set(&mut self, name: &str, set: Set) -> Result<()> { + let set = set.normalized().map_err(|e| context("set", name, e))?; + replace_existing(&mut self.sets, "set", name, set) + } + + pub fn remove_set(&mut self, name: &str) -> Result<()> { + not_in_use("set", name, self.set_users(name))?; + remove_existing(&mut self.sets, "set", name).map(|_| ()) + } + + /// The definitions that refer to a set, as `kind \`name\`` strings. + pub(crate) fn set_users(&self, set: &str) -> Vec { + let mut users = Vec::new(); + users.extend(self.sections.iter().filter(|(_, s)| s.set == set).map(|(n, _)| format!("section `{n}`"))); + users.extend(self.couplings.iter().filter(|(_, c)| c.set == set).map(|(n, _)| format!("coupling `{n}`"))); + users + .extend(self.boundary_conditions.iter().filter(|(_, b)| b.target().set() == Some(set)).map(|(n, _)| format!("boundary condition `{n}`"))); + users.extend(self.loads.iter().filter(|(_, l)| l.set() == Some(set)).map(|(n, _)| format!("load `{n}`"))); + users + } + + /// Resolve one set against the current mesh. The error names every member that is not on + /// the mesh. + pub fn resolve_set(&self, name: &str) -> Result { + let set = self.sets.get(name).ok_or_else(|| Error::NotFound(format!("set `{name}`")))?; + let (resolved, missing) = self.resolve_definition(set); + if missing.is_empty() { Ok(resolved) } else { Err(Error::NotFound(format!("set `{name}`: {}", missing.join("; ")))) } + } + + /// Resolve a definition, and list the members that are not on the mesh. + pub(crate) fn resolve_definition(&self, set: &Set) -> (ResolvedSet, Vec) { + let c = &self.classification; + let mut nodes = Vec::new(); + let mut elements = Vec::new(); + let mut faces = Vec::new(); + let mut missing = Vec::new(); + + // Geometry: mark the selected tags by role, then make one pass over the tagged elements. + let n = c.entities().len(); + let (mut regions, mut face_tags, mut edge_tags, mut vertex_tags) = (vec![false; n], vec![false; n], vec![false; n], vec![false; n]); + let mark = |marks: &mut Vec, tag: u32| { + if let Some(m) = marks.get_mut(tag as usize) { + *m = true; + } + }; + for r in &set.geometry { + let tags = c.tags_selected_by(r); + if tags.is_empty() { + missing.push(if c.is_empty() { + format!("{r} is not on the mesh (the mesh carries no geometry tags)") + } else { + format!("{r} is not on the mesh") + }); + continue; + } + for &t in &tags { + match r.entity.dimension() { + 0 => mark(&mut vertex_tags, t), + 1 => { + mark(&mut regions, t); + mark(&mut edge_tags, t); + } + 2 => { + mark(&mut regions, t); + mark(&mut face_tags, t); + } + _ => mark(&mut regions, t), + } + } + if r.entity.dimension() == 3 { + // A body's skin: every tagged face of the body. + for (i, e) in c.entities().iter().enumerate() { + if e.body == r.body && e.entity.dimension() == 2 { + mark(&mut face_tags, u32::try_from(i).unwrap_or(u32::MAX)); + } + } + } + } + let on = |marks: &[bool], tag: u32| marks.get(tag as usize).copied().unwrap_or(false); + if !set.geometry.is_empty() { + for (&label, tags) in c.elements() { + let Some(element) = self.elements.get(&label) else { continue }; + if tags.region.is_some_and(|t| on(®ions, t)) { + elements.push(label); + nodes.extend_from_slice(&element.nodes); + } + for (&face, &t) in &tags.faces { + if on(&face_tags, t) { + faces.push(ElementFace { element: label, face }); + } + } + for (&edge, &t) in &tags.edges { + if on(&edge_tags, t) { + let local = usize::from(edge).checked_sub(1).and_then(|i| element.kind.edges().get(i)).copied().unwrap_or_default(); + nodes.extend(local.iter().filter_map(|&i| element.nodes.get(i))); + } + } + } + nodes.extend(c.nodes().iter().filter(|(_, t)| on(&vertex_tags, **t)).map(|(n, _)| n)); + } + + // Explicit labels. + for &n in &set.nodes { + if self.nodes.contains_key(&n) { + nodes.push(n); + } else { + missing.push(format!("node {n} is not in the mesh")); + } + } + for &l in &set.elements { + match self.elements.get(&l) { + Some(e) => { + elements.push(l); + nodes.extend_from_slice(&e.nodes); + } + None => missing.push(format!("element {l} is not in the mesh")), + } + } + for &f in &set.faces { + match self.elements.get(&f.element) { + Some(e) if usize::from(f.face) <= e.kind.faces().len() => faces.push(f), + Some(e) => missing.push(format!("element {} ({}) has no face {}", f.element, e.kind.name(), f.face)), + None => missing.push(format!("element {} is not in the mesh", f.element)), + } + } + elements.sort_unstable(); + elements.dedup(); + faces.sort_unstable(); + faces.dedup(); + for f in &faces { + if let Some(e) = self.elements.get(&f.element) { + let local = usize::from(f.face).checked_sub(1).and_then(|i| e.kind.faces().get(i)).copied().unwrap_or_default(); + nodes.extend(local.iter().filter_map(|&i| e.nodes.get(i))); + } + } + nodes.sort_unstable(); + nodes.dedup(); + // Long lists of missing labels are summarised. + if missing.len() > 20 { + let more = missing.len() - 20; + missing.truncate(20); + missing.push(format!("and {more} more")); + } + (ResolvedSet { nodes, elements, faces }, missing) + } + + /// The node labels of an element face, in the face's order ([`crate::ElementKind::faces`]). + pub fn face_nodes(&self, face: ElementFace) -> Option> { + let e = self.elements.get(&face.element)?; + let local = e.kind.faces().get(usize::from(face.face).checked_sub(1)?)?; + local.iter().map(|&i| e.nodes.get(i).copied()).collect() + } +} diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs new file mode 100644 index 0000000..52fb47d --- /dev/null +++ b/crates/model/src/steps.rs @@ -0,0 +1,233 @@ +//! Analysis steps and the step × load/boundary-condition table. +//! +//! Steps run in order. Each step lists the boundary conditions and loads active in it, by name; +//! nothing carries over from one step to the next unless it is listed again, so the table a user +//! sees is the whole truth (no propagation, no hidden modified or inactive states). The state +//! before the first step is the unloaded, unconstrained model; a boundary condition active in +//! every step can be written at model level by formats that have one. + +use std::collections::BTreeSet; + +use serde::{Deserialize, Serialize}; + +use crate::{Error, Model, Result, check_name, context, finite, require}; + +/// Most steps in a model. +pub const MAX_STEPS: usize = 1000; + +/// Most modes a frequency step extracts. +pub const MAX_MODES: u32 = 100_000; + +/// What a step solves. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] +pub enum Procedure { + /// Linear static equilibrium under the active loads and boundary conditions. + Static {}, + /// Natural frequencies and mode shapes: the lowest `modes` modes, optionally only those + /// between `min_hz` and `max_hz` (cycles per time unit of the unit system). + Frequency { + modes: u32, + #[serde(default, skip_serializing_if = "Option::is_none")] + min_hz: Option, + #[serde(default, skip_serializing_if = "Option::is_none")] + max_hz: Option, + }, +} + +impl Procedure { + /// `static` or `frequency`. + pub fn name(&self) -> &'static str { + match self { + Procedure::Static {} => "static", + Procedure::Frequency { .. } => "frequency", + } + } + + pub fn validate(&self) -> Result<()> { + if let Procedure::Frequency { modes, min_hz, max_hz } = *self { + if !(1..=MAX_MODES).contains(&modes) { + return Err(Error::Invalid(format!("a frequency step extracts 1 to {MAX_MODES} modes, got {modes}"))); + } + for f in [min_hz, max_hz].into_iter().flatten() { + if finite("frequency", f)? < 0.0 { + return Err(Error::Invalid(format!("frequencies are not negative, got {f}"))); + } + } + if let (Some(lo), Some(hi)) = (min_hz, max_hz) + && lo >= hi + { + return Err(Error::Invalid(format!("the frequency range {lo}..{hi} is empty"))); + } + } + Ok(()) + } + + /// The output a step of this procedure requests unless told otherwise. + pub fn default_output(&self) -> BTreeSet { + match self { + Procedure::Static {} => OutputVariable::ALL.into_iter().collect(), + Procedure::Frequency { .. } => [OutputVariable::Displacement].into_iter().collect(), + } + } +} + +/// A result a step writes. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum OutputVariable { + /// Nodal displacements (mode shapes in a frequency step). + Displacement, + /// Reaction forces at constrained nodes. + ReactionForce, + Stress, + Strain, +} + +impl OutputVariable { + pub const ALL: [OutputVariable; 4] = + [OutputVariable::Displacement, OutputVariable::ReactionForce, OutputVariable::Stress, OutputVariable::Strain]; +} + +/// An analysis step. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Step { + pub name: String, + pub procedure: Procedure, + pub output: BTreeSet, + /// Boundary conditions active in this step. + #[serde(default)] + pub boundary_conditions: BTreeSet, + /// Loads active in this step. + #[serde(default)] + pub loads: BTreeSet, +} + +impl Step { + /// A step with the procedure's default output and nothing active. + pub fn new(name: impl Into, procedure: Procedure) -> Step { + Step { name: name.into(), output: procedure.default_output(), procedure, boundary_conditions: BTreeSet::new(), loads: BTreeSet::new() } + } +} + +impl Model { + /// The steps in analysis order. + pub fn steps(&self) -> &[Step] { + &self.steps + } + + pub fn step(&self, name: &str) -> Option<&Step> { + self.steps.iter().find(|s| s.name == name) + } + + fn step_index(&self, name: &str) -> Result { + self.steps.iter().position(|s| s.name == name).ok_or_else(|| Error::NotFound(format!("step `{name}`"))) + } + + fn step_mut(&mut self, name: &str) -> Result<&mut Step> { + self.steps.iter_mut().find(|s| s.name == name).ok_or_else(|| Error::NotFound(format!("step `{name}`"))) + } + + fn check_step(&self, step: &Step) -> Result<()> { + step.procedure.validate().map_err(|e| context("step", &step.name, e))?; + for b in &step.boundary_conditions { + require(&self.boundary_conditions, "boundary condition", b).map_err(|e| context("step", &step.name, e))?; + } + for l in &step.loads { + require(&self.loads, "load", l).map_err(|e| context("step", &step.name, e))?; + } + Ok(()) + } + + /// Add a step at the end, or just before the step named `before`. + pub fn add_step(&mut self, step: Step, before: Option<&str>) -> Result<()> { + check_name("step", &step.name)?; + if self.step(&step.name).is_some() { + return Err(Error::Invalid(format!("step `{}` already exists", step.name))); + } + if self.steps.len() >= MAX_STEPS { + return Err(Error::Limit(format!("a model holds at most {MAX_STEPS} steps"))); + } + self.check_step(&step)?; + let at = match before { + Some(b) => self.step_index(b)?, + None => self.steps.len(), + }; + self.steps.insert(at, step); + Ok(()) + } + + /// Replace a step's definition (procedure, output and table entries), keeping its place. + pub fn replace_step(&mut self, step: Step) -> Result<()> { + self.check_step(&step)?; + let name = step.name.clone(); + *self.step_mut(&name)? = step; + Ok(()) + } + + pub fn remove_step(&mut self, name: &str) -> Result<()> { + let i = self.step_index(name)?; + self.steps.remove(i); + Ok(()) + } + + /// Make a boundary condition active or inactive in a step. + pub fn set_boundary_condition_active(&mut self, step: &str, name: &str, active: bool) -> Result<()> { + require(&self.boundary_conditions, "boundary condition", name)?; + switch(&mut self.step_mut(step)?.boundary_conditions, name, active); + Ok(()) + } + + /// Make a load active or inactive in a step. + pub fn set_load_active(&mut self, step: &str, name: &str, active: bool) -> Result<()> { + require(&self.loads, "load", name)?; + switch(&mut self.step_mut(step)?.loads, name, active); + Ok(()) + } +} + +fn switch(active: &mut BTreeSet, name: &str, on: bool) { + if on { + active.insert(name.to_string()); + } else { + active.remove(name); + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{BoundaryCondition, Set, Target}; + + #[test] + fn steps_keep_order_and_an_explicit_table() { + let mut m = Model::default(); + m.add_set("base", Set { nodes: vec![1], ..Set::default() }).unwrap(); + m.add_boundary_condition("fix", BoundaryCondition::Fixed { target: Target::Set("base".into()) }).unwrap(); + m.add_step(Step::new("Load", Procedure::Static {}), None).unwrap(); + m.add_step(Step::new("Modes", Procedure::Frequency { modes: 10, min_hz: None, max_hz: None }), Some("Load")).unwrap(); + assert_eq!(m.steps().iter().map(|s| s.name.as_str()).collect::>(), ["Modes", "Load"]); + assert!(m.add_step(Step::new("Load", Procedure::Static {}), None).is_err()); + assert!(m.add_step(Step::new("X", Procedure::Static {}), Some("Nope")).is_err()); + assert!(m.add_step(Step::new("Y", Procedure::Frequency { modes: 0, min_hz: None, max_hz: None }), None).is_err()); + assert!(m.add_step(Step::new("Y", Procedure::Frequency { modes: 5, min_hz: Some(10.0), max_hz: Some(1.0) }), None).is_err()); + assert!(serde_json::from_value::(serde_json::json!({"type": "static", "modes": 3})).is_err(), "typos are errors"); + + m.set_boundary_condition_active("Load", "fix", true).unwrap(); + assert!(m.step("Load").unwrap().boundary_conditions.contains("fix")); + assert!(m.step("Modes").unwrap().boundary_conditions.is_empty(), "no propagation"); + assert!(m.set_boundary_condition_active("Load", "nope", true).is_err()); + assert!(m.set_load_active("Load", "fix", true).is_err(), "loads and boundary conditions are separate"); + + // Nothing still referred to can be deleted: the set by the boundary condition, the + // boundary condition by the step that lists it. + assert!(matches!(m.remove_set("base"), Err(Error::InUse(_)))); + assert!(m.remove_boundary_condition("fix").unwrap_err().to_string().contains("step `Load`")); + m.set_boundary_condition_active("Load", "fix", false).unwrap(); + m.remove_boundary_condition("fix").unwrap(); + m.remove_set("base").unwrap(); + m.remove_step("Modes").unwrap(); + assert_eq!(m.steps().len(), 1); + } +} diff --git a/crates/model/src/units.rs b/crates/model/src/units.rs new file mode 100644 index 0000000..f75fe47 --- /dev/null +++ b/crates/model/src/units.rs @@ -0,0 +1,146 @@ +//! The unit system label. +//! +//! The model never converts: numbers are stored as entered. The label records which consistent +//! system they are in, so the check can warn when a magnitude is implausible for it (a steel +//! modulus of 210 000 in a metre–pascal model, a density of 7850 in a millimetre–tonne one), and +//! so writers can say which system a deck is in. + +use serde::{Deserialize, Serialize}; + +/// A consistent unit system: length, mass, time, force, stress. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +pub enum UnitSystem { + /// Not chosen: magnitudes are not checked. + #[default] + #[serde(rename = "unspecified")] + Unspecified, + /// Metre, kilogram, second, newton, pascal. + #[serde(rename = "m-kg-s-N-Pa")] + MKgS, + /// Millimetre, tonne, second, newton, megapascal: the usual choice for parts in millimetres. + #[serde(rename = "mm-t-s-N-MPa")] + MmTonneS, + /// Millimetre, kilogram, millisecond, kilonewton, gigapascal. + #[serde(rename = "mm-kg-ms-kN-GPa")] + MmKgMs, + /// Millimetre, gram, millisecond, newton, megapascal. + #[serde(rename = "mm-g-ms-N-MPa")] + MmGMs, + /// Inch, pound-force second squared per inch, second, pound-force, psi. + #[serde(rename = "in-lbf-s-psi")] + InLbfS, +} + +impl UnitSystem { + pub const ALL: [UnitSystem; 6] = + [UnitSystem::Unspecified, UnitSystem::MKgS, UnitSystem::MmTonneS, UnitSystem::MmKgMs, UnitSystem::MmGMs, UnitSystem::InLbfS]; + + /// The label used in commands and files (`mm-t-s-N-MPa`, …). + pub fn label(self) -> &'static str { + match self { + UnitSystem::Unspecified => "unspecified", + UnitSystem::MKgS => "m-kg-s-N-Pa", + UnitSystem::MmTonneS => "mm-t-s-N-MPa", + UnitSystem::MmKgMs => "mm-kg-ms-kN-GPa", + UnitSystem::MmGMs => "mm-g-ms-N-MPa", + UnitSystem::InLbfS => "in-lbf-s-psi", + } + } + + pub fn from_label(label: &str) -> Option { + UnitSystem::ALL.into_iter().find(|u| u.label() == label) + } + + /// The length unit (`m`, `mm`, `in`), or `None` when unspecified. + pub fn length(self) -> Option<&'static str> { + match self { + UnitSystem::Unspecified => None, + UnitSystem::MKgS => Some("m"), + UnitSystem::MmTonneS | UnitSystem::MmKgMs | UnitSystem::MmGMs => Some("mm"), + UnitSystem::InLbfS => Some("in"), + } + } + + /// Size of one stress unit in pascals, one density unit in kg/m³ and one acceleration unit in + /// m/s², or `None` when unspecified. + fn scales(self) -> Option { + const IN: f64 = 0.0254; + const LBF: f64 = 4.448_221_615_260_5; + Some(match self { + UnitSystem::Unspecified => return None, + UnitSystem::MKgS => Scales { stress: 1.0, density: 1.0, acceleration: 1.0 }, + UnitSystem::MmTonneS => Scales { stress: 1e6, density: 1e12, acceleration: 1e-3 }, + UnitSystem::MmKgMs => Scales { stress: 1e9, density: 1e9, acceleration: 1e3 }, + UnitSystem::MmGMs => Scales { stress: 1e6, density: 1e6, acceleration: 1e3 }, + // Mass unit lbf·s²/in = LBF / IN kg; density unit = that per cubic inch. + UnitSystem::InLbfS => Scales { stress: LBF / (IN * IN), density: LBF / IN / (IN * IN * IN), acceleration: IN }, + }) + } + + /// A warning when a Young's or shear modulus is implausible for any engineering material in + /// this system (outside 10 kPa – 2 TPa). + pub fn check_modulus(self, what: &str, value: f64) -> Option { + self.check(what, value, |s| s.stress, (1e4, 2e12), "Pa") + } + + /// A warning when a mass density is implausible (outside 1 – 25 000 kg/m³). + pub fn check_density(self, value: f64) -> Option { + self.check("density", value, |s| s.density, (1.0, 25_000.0), "kg/m³") + } + + /// A warning when a gravitational acceleration is implausible (outside 0.1 – 100 m/s²). + pub fn check_gravity(self, magnitude: f64) -> Option { + self.check("gravity", magnitude, |s| s.acceleration, (0.1, 100.0), "m/s²") + } + + fn check(self, what: &str, value: f64, scale: fn(&Scales) -> f64, (lo, hi): (f64, f64), si: &str) -> Option { + let s = self.scales()?; + let si_value = value.abs() * scale(&s); + if (lo..=hi).contains(&si_value) { + return None; + } + let k = scale(&s); + Some(format!( + "{what} {value} is {si_value:.3e} {si}, outside the plausible {lo:.0e}–{hi:.0e} {si} ({:.3e}–{:.3e} in {}); check the unit system", + lo / k, + hi / k, + self.label() + )) + } +} + +struct Scales { + stress: f64, + density: f64, + acceleration: f64, +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn steel_is_plausible_in_every_system_and_wrong_systems_warn() { + // Steel: E = 210 GPa, ρ = 7850 kg/m³, g = 9.81 m/s², in each system's own units. + let steel = [ + (UnitSystem::MKgS, 2.1e11, 7850.0, 9.81), + (UnitSystem::MmTonneS, 2.1e5, 7.85e-9, 9810.0), + (UnitSystem::MmKgMs, 210.0, 7.85e-6, 9.81e-3), + (UnitSystem::MmGMs, 2.1e5, 7.85e-3, 9.81e-3), + (UnitSystem::InLbfS, 3.05e7, 7.35e-4, 386.1), + ]; + for (u, e, rho, g) in steel { + assert_eq!(u.check_modulus("E", e), None, "{u:?}"); + assert_eq!(u.check_density(rho), None, "{u:?}"); + assert_eq!(u.check_gravity(g), None, "{u:?}"); + assert_eq!(UnitSystem::from_label(u.label()), Some(u)); + assert_eq!(serde_json::to_value(u).unwrap(), serde_json::json!(u.label())); + } + // SI numbers in a millimetre–tonne model. + assert!(UnitSystem::MmTonneS.check_modulus("E", 2.1e11).is_some()); + assert!(UnitSystem::MmTonneS.check_density(7850.0).is_some()); + assert!(UnitSystem::MKgS.check_gravity(9810.0).is_some()); + assert_eq!(UnitSystem::Unspecified.check_density(1e30), None); + assert!(UnitSystem::MKgS.check_density(f64::NAN).is_some()); + } +} diff --git a/crates/model/tests/definition.rs b/crates/model/tests/definition.rs new file mode 100644 index 0000000..1032105 --- /dev/null +++ b/crates/model/tests/definition.rs @@ -0,0 +1,396 @@ +//! A whole model definition on a tagged mesh: sets by geometry name survive a remesh, the check +//! names what is missing, and the model round-trips through JSON. + +use modelcraft_model::*; + +/// A bar of `n` unit hexahedra along X, labels starting after `offset`, tagged as body `Bar` with +/// faces `left` (x = 0) and `right` (x = n), and the vertex `tip` at (n, 0, 0). +fn bar(m: &mut Model, n: u32, offset: u32) { + m.clear_geometry_tags(); + for e in m.elements().keys().copied().collect::>() { + m.remove_element(e).unwrap(); + } + for v in m.nodes().keys().copied().collect::>() { + m.remove_node(v).unwrap(); + } + let node = |i: u32, j: u32, k: u32| offset + 1 + i * 4 + j * 2 + k; + for i in 0..=n { + for j in 0..2 { + for k in 0..2 { + m.set_node(node(i, j, k), [f64::from(i), f64::from(j), f64::from(k)]).unwrap(); + } + } + } + let body = m.geometry_tag(&GeometryRef::body("Bar").unwrap()).unwrap(); + let left = m.geometry_tag(&GeometryRef::face("Bar", "left").unwrap()).unwrap(); + let right = m.geometry_tag(&GeometryRef::face("Bar", "right").unwrap()).unwrap(); + let tip = m.geometry_tag(&GeometryRef::new("Bar", Entity::Vertex("tip".into())).unwrap()).unwrap(); + for i in 0..n { + let label = offset + 1 + i; + let nodes = vec![ + node(i, 0, 0), + node(i + 1, 0, 0), + node(i + 1, 1, 0), + node(i, 1, 0), + node(i, 0, 1), + node(i + 1, 0, 1), + node(i + 1, 1, 1), + node(i, 1, 1), + ]; + m.set_element(label, ElementKind::Hex8, nodes).unwrap(); + m.tag_element(label, body).unwrap(); + if i == 0 { + m.tag_element_face(label, 6, left).unwrap(); + } + if i + 1 == n { + m.tag_element_face(label, 4, right).unwrap(); + } + } + m.tag_node(node(n, 0, 0), tip).unwrap(); +} + +fn defined() -> Model { + let mut m = Model::default(); + bar(&mut m, 2, 0); + m.set_units(UnitSystem::MmTonneS); + m.add_set("all", Set::of_geometry(vec![GeometryRef::body("Bar").unwrap()])).unwrap(); + m.add_set("root", Set::of_geometry(vec![GeometryRef::face("Bar", "left").unwrap()])).unwrap(); + m.add_set("end", Set::of_geometry(vec![GeometryRef::face("Bar", "right").unwrap()])).unwrap(); + m.add_material("Steel", Material { elastic: Some(Elastic::Isotropic { young: 2.1e5, poisson: 0.3 }), density: Some(7.85e-9) }).unwrap(); + m.add_section("bar", Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() }) + .unwrap(); + m.add_reference_point("RP", ReferencePoint { xyz: [3.0, 0.5, 0.5], label: 0 }).unwrap(); + m.add_coupling("load-pad", Coupling { reference_point: "RP".into(), set: "end".into(), kind: CouplingKind::Distributing, dofs: Dofs::all() }) + .unwrap(); + m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("root".into()) }).unwrap(); + m.add_load("push", Load::Force { target: Target::ReferencePoint("RP".into()), vector: [0.0, 0.0, -100.0], distribution: Distribution::Total }) + .unwrap(); + m.add_load("press", Load::Pressure { set: "end".into(), magnitude: 1.0 }).unwrap(); + let mut s = Step::new("Bend", Procedure::Static {}); + s.boundary_conditions.insert("clamp".into()); + s.loads.extend(["push".to_string(), "press".to_string()]); + m.add_step(s, None).unwrap(); + let mut f = Step::new("Modes", Procedure::Frequency { modes: 6, min_hz: None, max_hz: None }); + f.boundary_conditions.insert("clamp".into()); + m.add_step(f, None).unwrap(); + m +} + +#[test] +fn a_complete_model_resolves_and_survives_a_remesh() { + let mut m = defined(); + assert_eq!(m.check(), vec![]); + let r = m.resolve().unwrap(); + assert_eq!(r.sets["all"].elements, vec![1, 2]); + assert_eq!(r.sets["all"].nodes.len(), 12); + assert_eq!(r.sets["all"].faces, vec![ElementFace { element: 1, face: 6 }, ElementFace { element: 2, face: 4 }], "a body's skin"); + assert_eq!(r.sets["root"].faces, vec![ElementFace { element: 1, face: 6 }]); + assert_eq!(r.sets["root"].nodes, vec![1, 2, 3, 4]); + assert!(r.sets["root"].elements.is_empty(), "a face holds no solid elements"); + assert_eq!(r.sets["end"].nodes, vec![9, 10, 11, 12]); + assert_eq!(r.reference_point_nodes["RP"], MAX_LABEL, "from the top of the label range"); + assert_eq!(r.nodal_loads["push"], vec![(MAX_LABEL, [0.0, 0.0, -100.0])]); + assert_eq!(r.target_nodes(&Target::ReferencePoint("RP".into())), [MAX_LABEL]); + + // The pressure acts on the end face's outward-oriented nodes (x = 2). + let face = m.face_nodes(ElementFace { element: 2, face: 4 }).unwrap(); + assert!(face.iter().all(|n| m.node(*n).unwrap().xyz[0] == 2.0)); + + // Remesh with four elements and other labels: the same definitions resolve on the new mesh. + bar(&mut m, 4, 100); + assert_eq!(m.check(), vec![]); + let r = m.resolve().unwrap(); + assert_eq!(r.sets["all"].elements, vec![101, 102, 103, 104]); + assert_eq!(r.sets["end"].faces, vec![ElementFace { element: 104, face: 4 }]); + assert_eq!(r.sets["end"].nodes, vec![117, 118, 119, 120]); + assert_eq!(r.reference_point_nodes["RP"], MAX_LABEL, "a remesh does not renumber reference points"); + + // Vertex members resolve to their node. + m.add_set("tip", Set::of_geometry(vec![GeometryRef::new("Bar", Entity::Vertex("tip".into())).unwrap()])).unwrap(); + assert_eq!(m.resolve_set("tip").unwrap().nodes, vec![117]); +} + +#[test] +fn the_check_names_what_is_missing() { + let mut m = defined(); + // A set naming a face the mesh does not carry, and explicit labels not in the mesh. + m.add_set("ghost", Set { geometry: vec![GeometryRef::face("Bar", "top").unwrap()], nodes: vec![999], ..Set::default() }).unwrap(); + let e = m.resolve_set("ghost").unwrap_err().to_string(); + assert!(e.contains("face top of body Bar is not on the mesh") && e.contains("node 999"), "{e}"); + m.add_set("side", Set { faces: vec![ElementFace { element: 1, face: 3 }], ..Set::default() }).unwrap(); + m.add_load("ghost-load", Load::Pressure { set: "ghost".into(), magnitude: 1.0 }).unwrap(); + m.add_load("twist", Load::Moment { target: Target::Set("end".into()), vector: [1.0, 0.0, 0.0], distribution: Distribution::Total }).unwrap(); + m.add_material("Soft", Material { elastic: Some(Elastic::Isotropic { young: 2.1e11, poisson: 0.3 }), density: None }).unwrap(); + m.add_section("again", Section { set: "side".into(), material: "Soft".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() }) + .unwrap(); + m.add_section( + "shell", + Section { + set: "root".into(), + material: "Steel".into(), + kind: SectionKind::Shell { thickness: 1.0, offset: 0.0, points: 5 }, + formulation: Formulation::default(), + }, + ) + .unwrap(); + + let diagnostics = m.check(); + let text: Vec = diagnostics.iter().map(ToString::to_string).collect(); + let has = |needle: &str| text.iter().any(|t| t.contains(needle)); + assert!(has("error: set `ghost`: face top of body Bar is not on the mesh"), "{text:#?}"); + assert!(has("error: set `ghost`: node 999 is not in the mesh"), "{text:#?}"); + assert!(has("error: load `ghost-load`: set `ghost` has no element faces"), "{text:#?}"); + assert!(has("error: load `twist`: set `end` has no nodes with rotations"), "{text:#?}"); + assert!(has("error: section `shell`: set `root` has no elements"), "{text:#?}"); + assert!(has("error: section `again`: set `side` has no elements"), "{text:#?}"); + assert!(has("error: material `Soft`: has no density, and step `Modes` extracts frequencies"), "{text:#?}"); + assert!(has("warning: material `Soft`: Young's modulus 210000000000"), "{text:#?}"); + assert!(has("warning: load `twist`: is not active in any step"), "{text:#?}"); + assert_eq!(diagnostics.first().map(|d| d.severity), Some(Severity::Error), "errors first"); + assert!(m.resolve().is_err()); + + // Two sections over the same elements, and an element with none. + let mut m = defined(); + m.add_section("twice", Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() }) + .unwrap(); + m.set_node(50, [9.0, 9.0, 9.0]).unwrap(); + m.set_node(51, [10.0, 9.0, 9.0]).unwrap(); + m.set_element(9, ElementKind::Rod2, vec![50, 51]).unwrap(); + let text: Vec = m.check().iter().map(ToString::to_string).collect(); + assert!(text.iter().any(|t| t.contains("section `twice`: elements 1, 2 are also in section `bar`")), "{text:#?}"); + assert!(text.iter().any(|t| t.contains("model: elements 9 have no section")), "{text:#?}"); + + // Incompatible modes on tetrahedra, a beam oriented along its axis. + let mut m = defined(); + m.replace_section( + "bar", + Section { + set: "all".into(), + material: "Steel".into(), + kind: SectionKind::Solid {}, + formulation: Formulation { incompatible_modes: true, ..Formulation::default() }, + }, + ) + .unwrap(); + assert_eq!(m.check().iter().filter(|d| d.severity == Severity::Error).count(), 0, "hex8 takes incompatible modes"); + m.set_element(2, ElementKind::Tet4, vec![5, 6, 7, 9]).unwrap(); + let body = m.geometry_tag(&GeometryRef::body("Bar").unwrap()).unwrap(); + m.tag_element(2, body).unwrap(); + assert!(m.check().iter().any(|d| d.message.contains("incompatible modes apply to hex8 elements only, not elements 2"))); + m.set_element(30, ElementKind::Beam2, vec![1, 5]).unwrap(); + m.add_set("beam", Set { elements: vec![30], ..Set::default() }).unwrap(); + let profile = BeamProfile::Rectangle { width: 1.0, height: 2.0 }; + m.add_section( + "beam", + Section { + set: "beam".into(), + material: "Steel".into(), + kind: SectionKind::Beam { profile, orientation: [-3.0, 0.0, 0.0] }, + formulation: Formulation::default(), + }, + ) + .unwrap(); + assert!(m.check().iter().any(|d| d.subject == "section `beam`" && d.message.contains("parallel to elements 30"))); +} + +#[test] +fn definitions_refer_only_to_existing_items_and_cannot_be_deleted_while_used() { + let mut m = defined(); + let section = |set: &str, material: &str| Section { + set: set.into(), + material: material.into(), + kind: SectionKind::Solid {}, + formulation: Formulation::default(), + }; + assert!(matches!(m.add_section("x", section("nope", "Steel")), Err(Error::NotFound(_)))); + assert!(matches!(m.add_section("x", section("all", "Nope")), Err(Error::NotFound(_)))); + assert!( + m.add_load("x", Load::Force { target: Target::ReferencePoint("nope".into()), vector: [1.0, 0.0, 0.0], distribution: Distribution::Total }) + .is_err() + ); + assert!( + m.add_coupling("x", Coupling { reference_point: "RP".into(), set: "end".into(), kind: CouplingKind::Kinematic, dofs: Dofs::none() }).is_err() + ); + + let e = m.remove_set("end").unwrap_err().to_string(); + assert!(e.contains("coupling `load-pad`") && e.contains("load `press`"), "{e}"); + assert!(m.remove_material("Steel").unwrap_err().to_string().contains("section `bar`")); + assert!(m.remove_reference_point("RP").unwrap_err().to_string().contains("load `push`")); + assert!(m.remove_load("push").unwrap_err().to_string().contains("step `Bend`")); + m.set_load_active("Bend", "push", false).unwrap(); + m.remove_load("push").unwrap(); + m.remove_coupling("load-pad").unwrap(); + m.remove_reference_point("RP").unwrap(); + + // Removing a tagged element drops its tags; a new element under the label starts untagged. + m.remove_element(1).unwrap(); + assert!(!m.classification().elements().contains_key(&1)); +} + +#[test] +fn the_model_round_trips_through_json() { + let m = defined(); + let json = serde_json::to_string(&m).unwrap(); + let back: Model = serde_json::from_str(&json).unwrap(); + assert_eq!(back, m); + assert_eq!(serde_json::to_string(&back).unwrap(), json, "deterministic"); + // Reading checks what commands check: a dangling reference, a bad tag, a typo. + let mut v = serde_json::to_value(&m).unwrap(); + v["sections"]["bar"]["material"] = serde_json::json!("Nope"); + assert!(serde_json::from_value::(v).unwrap_err().to_string().contains("material `Nope`")); + let mut v = serde_json::to_value(&m).unwrap(); + v["classification"]["elements"]["1"]["faces"]["6"] = serde_json::json!(0); + assert!(serde_json::from_value::(v).is_err(), "tag 0 is the body, not a face"); + let mut v = serde_json::to_value(&m).unwrap(); + v["materials"]["Steel"]["elastic"]["young"] = serde_json::json!(-1.0); + assert!(serde_json::from_value::(v).is_err()); + let mut v = serde_json::to_value(&m).unwrap(); + v["stepz"] = serde_json::json!([]); + assert!(serde_json::from_value::(v).is_err()); + + // A model file from before the definitions existed still reads. + let old: Model = serde_json::from_value(serde_json::json!({"nodes": {}, "elements": {}})).unwrap(); + assert_eq!(old, Model::default()); +} + +#[test] +fn tags_are_checked() { + let mut m = defined(); + let face = m.geometry_tag(&GeometryRef::face("Bar", "left").unwrap()).unwrap(); + let body = m.geometry_tag(&GeometryRef::body("Bar").unwrap()).unwrap(); + assert!(m.tag_element(1, face).is_err(), "a solid fills a body, not a face"); + assert!(m.tag_element_face(1, 7, face).is_err()); + assert!(m.tag_element_face(1, 0, face).is_err()); + assert!(m.tag_element_face(1, 1, body).is_err()); + assert!(m.tag_element_edge(1, 13, 0).is_err()); + assert!(m.tag_node(999, 0).is_err()); + assert!(m.tag_element(1, 9999).is_err()); + assert!(m.tag_element(77, body).is_err()); + assert_eq!(m.geometry_tag(&GeometryRef::face("Bar", "left").unwrap()).unwrap(), face, "interned once"); +} + +/// A block `[0, 1] × [0, 1] × [0, 1]` of `n × n × n` hex8 with graded spacing (so area weighting +/// matters), tagged as body `Block` with its face `end` at x = 1. +fn block(m: &mut Model, n: u32) { + m.clear_geometry_tags(); + for e in m.elements().keys().copied().collect::>() { + m.remove_element(e).unwrap(); + } + for v in m.nodes().keys().copied().collect::>() { + m.remove_node(v).unwrap(); + } + let at = |i: u32| (f64::from(i) / f64::from(n)).powi(2); + let node = |i: u32, j: u32, k: u32| 1 + i + j * (n + 1) + k * (n + 1) * (n + 1); + for k in 0..=n { + for j in 0..=n { + for i in 0..=n { + m.set_node(node(i, j, k), [at(i), at(j), at(k)]).unwrap(); + } + } + } + let body = m.geometry_tag(&GeometryRef::body("Block").unwrap()).unwrap(); + let end = m.geometry_tag(&GeometryRef::face("Block", "end").unwrap()).unwrap(); + let mut label = 0; + for k in 0..n { + for j in 0..n { + for i in 0..n { + label += 1; + let nodes = vec![ + node(i, j, k), + node(i + 1, j, k), + node(i + 1, j + 1, k), + node(i, j + 1, k), + node(i, j, k + 1), + node(i + 1, j, k + 1), + node(i + 1, j + 1, k + 1), + node(i, j + 1, k + 1), + ]; + m.set_element(label, ElementKind::Hex8, nodes).unwrap(); + m.tag_element(label, body).unwrap(); + if i + 1 == n { + m.tag_element_face(label, 4, end).unwrap(); + } + } + } + } +} + +fn sum(values: &[(u32, [f64; 3])]) -> [f64; 3] { + values.iter().fold([0.0; 3], |s, (_, v)| [s[0] + v[0], s[1] + v[1], s[2] + v[2]]) +} + +#[test] +fn a_total_force_on_a_face_is_the_same_on_every_mesh() { + let mut m = Model::default(); + block(&mut m, 2); + m.set_units(UnitSystem::MmTonneS); + m.add_set("all", Set::of_geometry(vec![GeometryRef::body("Block").unwrap()])).unwrap(); + m.add_set("end", Set::of_geometry(vec![GeometryRef::face("Block", "end").unwrap()])).unwrap(); + m.add_material("Steel", Material { elastic: Some(Elastic::Isotropic { young: 2.1e5, poisson: 0.3 }), density: Some(7.85e-9) }).unwrap(); + m.add_section("block", Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() }) + .unwrap(); + let force = |distribution| Load::Force { target: Target::Set("end".into()), vector: [100.0, 0.0, -50.0], distribution }; + m.add_load("total", force(Distribution::Total)).unwrap(); + m.add_load("each", force(Distribution::PerNode)).unwrap(); + m.add_load("drag", Load::Traction { set: "end".into(), vector: [0.0, 2.0, 0.0] }).unwrap(); + let mut s = Step::new("Pull", Procedure::Static {}); + s.loads.extend(["total".to_string(), "each".to_string(), "drag".to_string()]); + m.add_step(s, None).unwrap(); + + let close = |a: [f64; 3], b: [f64; 3]| a.iter().zip(b).all(|(x, y)| (x - y).abs() < 1e-9); + let coarse = m.resolve().unwrap(); + block(&mut m, 5); + let fine = m.resolve().unwrap(); + for r in [&coarse, &fine] { + assert!(close(sum(&r.nodal_loads["total"]), [100.0, 0.0, -50.0]), "{:?}", sum(&r.nodal_loads["total"])); + // Traction × area (the end face is 1 × 1). + assert!(close(sum(&r.nodal_loads["drag"]), [0.0, 2.0, 0.0])); + assert_eq!(r.nodal_loads["total"].len(), r.sets["end"].nodes.len()); + } + // A per-node force multiplies with the node count, by request. + assert!(close(sum(&coarse.nodal_loads["each"]), [900.0, 0.0, -450.0])); + assert!(close(sum(&fine.nodal_loads["each"]), [3600.0, 0.0, -1800.0])); + let weights: Vec = fine.nodal_loads["total"].iter().map(|(_, v)| v[0]).collect(); + assert!(weights.iter().any(|w| (w - weights[0]).abs() > 1e-6), "graded faces give unequal shares"); + + // A node set without faces splits a total force equally. + m.add_set("pins", Set { nodes: vec![1, 2, 3, 4], ..Set::default() }).unwrap(); + m.add_load("pinned", Load::Force { target: Target::Set("pins".into()), vector: [4.0, 0.0, 0.0], distribution: Distribution::Total }).unwrap(); + assert_eq!( + m.resolve().unwrap().nodal_loads["pinned"], + vec![(1, [1.0, 0.0, 0.0]), (2, [1.0, 0.0, 0.0]), (3, [1.0, 0.0, 0.0]), (4, [1.0, 0.0, 0.0])] + ); + + // A traction needs faces. + m.add_load("nowhere", Load::Traction { set: "pins".into(), vector: [1.0, 0.0, 0.0] }).unwrap(); + assert!(m.check().iter().any(|d| d.subject == "load `nowhere`" && d.message.contains("no element faces"))); +} + +#[test] +fn reference_point_labels_are_stable() { + let mut m = defined(); + m.add_reference_point("A", ReferencePoint { xyz: [0.0; 3], label: 0 }).unwrap(); + m.add_reference_point("Z", ReferencePoint { xyz: [0.0; 3], label: 0 }).unwrap(); + let before = m.reference_point_nodes(); + assert_eq!((before["RP"], before["A"], before["Z"]), (MAX_LABEL, MAX_LABEL - 1, MAX_LABEL - 2), "in creation order"); + // Another point, more nodes, a remesh and an edit renumber nothing. + m.add_reference_point("B", ReferencePoint { xyz: [1.0; 3], label: 0 }).unwrap(); + m.set_node(5000, [9.0; 3]).unwrap(); + bar(&mut m, 3, 40); + m.replace_reference_point("A", ReferencePoint { xyz: [2.0; 3], label: 0 }).unwrap(); + let after = m.reference_point_nodes(); + for name in ["RP", "A", "Z"] { + assert_eq!(after[name], before[name], "{name}"); + } + // Removing one frees its label for the next; explicit labels are checked. + m.remove_reference_point("Z").unwrap(); + m.add_reference_point("C", ReferencePoint { xyz: [0.0; 3], label: 0 }).unwrap(); + assert_eq!(m.reference_point_nodes()["C"], MAX_LABEL - 2); + assert!(m.add_reference_point("D", ReferencePoint { xyz: [0.0; 3], label: MAX_LABEL }).is_err(), "taken"); + assert!(m.add_reference_point("D", ReferencePoint { xyz: [0.0; 3], label: MAX_LABEL + 1 }).is_err()); + m.add_reference_point("D", ReferencePoint { xyz: [0.0; 3], label: 41 }).unwrap(); + assert!( + m.check().iter().any(|d| d.subject == "reference point `D`" && d.message.contains("also a mesh node")), + "41 is a node of the remeshed bar" + ); +} diff --git a/docs/architecture.md b/docs/architecture.md index f231cf7..d145885 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -18,7 +18,7 @@ crates in strictly lower layers. | Layer | Crate | Role | Status | |---|---|---|---| -| L0 | `modelcraft-model` | Solver-neutral FE model: nodes, elements by label; materials, sections, sets, BCs, loads and steps as they land. Knows no file format. | started | +| L0 | `modelcraft-model` | Solver-neutral FE model: nodes and elements by label, their geometry tags, sets, materials, sections, couplings, BCs, loads and steps, the model check and the resolved view the writers read. Knows no file format and no geometry type. | started | | L1 | `modelcraft-geometry` | The only door to SolidCraft: opens `.solidcraft` designs (rebuilt from their timeline, bodies with persistent face names), B-rep bodies, measures, tessellation into a face-tagged `SurfaceMesh`; face-tagged 3MF files read into the same type without the kernel (docs/formats/tagged-mesh.md). No SolidCraft type in its API. | started | | L1 | `modelcraft-mesh` | Surface preparation (weld, watertightness, feature edges), robust predicates, Delaunay tetrahedralization and refinement, Voronoi-based optimization, hex meshing, quality metrics. | started | | L1 | `modelcraft-solver` | Element formulations, assembly, sparse linear static and modal solution. | planned | @@ -78,6 +78,48 @@ Rules beyond the layer order: 6. **Quality.** Scaled Jacobian, aspect ratio, dihedral angles; every generated mesh is checked for positive Jacobians and conserved volume before it reaches the model. +## Model definition + +`modelcraft-model` holds the analysis model in no solver's vocabulary (the writers map it to +element names, cards and keywords): + +- **Geometry tags.** Geometry is referred to by name only: a body and, inside it, the persistent + name of a face, edge or vertex (`{"body": "Body1", "face": "F2:top"}`). A mesher tags what it + makes: each element with the body (solids), face (shells) or edge (beams) it fills, each + boundary element face with its geometry face, element edges with geometry edges, nodes with + vertices (`Model::geometry_tag`, `tag_element`, `tag_element_face`, …; `mesh.tag` by hand). + Replacing or removing an element drops its tags; a remesh clears them and writes new ones. +- **Sets** are defined by geometry names, explicit labels or both, and resolve against the + current mesh into sorted nodes, elements and element faces. A face name also selects the pieces + a kernel splits it into (`name#0`, `name#1`). Because the definition names geometry, a remesh, + or an edit of the design followed by a remesh, keeps every set, section and load in place. +- **Element faces** are numbered per kind (`ElementKind::faces`): solids' faces are oriented with + the right-hand normal pointing out; shells have a positive and a negative side. +- **Definitions** (materials, sections with their element formulation, reference points, + couplings, boundary conditions, loads, steps) refer to each other by name. A reference must + exist when it is made, and a definition still referred to (a load listed by a step included) + cannot be deleted; what depends on + the mesh is checked when the model is resolved. +- **Steps** list the boundary conditions and loads active in them. Nothing propagates from one + step to the next. +- **Units** are a label (`mm-t-s-N-MPa`, …); numbers are never converted, and the check warns + when a modulus, density or gravity is implausible for the chosen system. +- **Check and resolve.** `Model::check` lists errors and warnings, each naming its item (a set + member not on the mesh, an element with no section or in two, a moment on solid nodes, a + missing density for a frequency step, …). `Model::resolve` returns what writers and the solver + read besides the definitions: every set resolved (each element then lies in exactly one + section's set), the node label of every reference point, and every force, moment and traction + as nodal values. It fails when the check finds an error. +- **Loads that keep their meaning on any mesh.** A force or moment on a set is a total by + default: split by the consistent nodal weights `∫ N_i dA` of the set's faces (Gauss quadrature + on each face's own linear or quadratic geometry), or equally over a set of nodes only. Applying + it at every node is an explicit, mesh-dependent opt-in (`per_node`). A traction is force per + area on faces. The expansion lives in `resolve`, so the deck writers and the solver share it. +- **Reference point labels** are stored with the point: chosen from the top of the label range + when it is created, so remeshing or adding nodes and points never renumbers them. +- **Files.** A model deserializes by replaying its contents through the same checked methods the + commands use, so a file cannot hold what commands would refuse. + ## Commands Every user-visible action is a `CommandSpec` in `crates/engine/src/cmd/.rs` with a JSON diff --git a/docs/parity.md b/docs/parity.md index a74b65b..baf5de4 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -42,16 +42,16 @@ products are run, no documentation text is copied). Status: **done**, **started* | Mesh | Hex-dominant by recombination; polycube | M6 | yes | planned | | Mesh | Smoothing, untangling, boundary padding, layers through thin walls | M6 | yes | planned | | Mesh | Mesh editing (move, merge, split, collapse, remesh region), imported orphan meshes | M6 | yes | started (node/element add/remove) | -| Model | Sets and surfaces bound to geometry names | M3 | yes | planned | -| Model | Materials: elastic (isotropic, orthotropic), density, unit-system check | M3 | | planned | +| Model | Sets and surfaces bound to geometry names | M3 | yes | done | +| Model | Materials: elastic (isotropic, orthotropic), density, unit-system check | M3 | | done | | Model | Materials: plasticity, hyperelasticity, damage, thermal | M8 | | planned | -| Model | Sections: solid, shell, beam profiles; assignment by set | M3 | | planned | +| Model | Sections: solid, shell, beam profiles; assignment by set | M3 | | done | | Model | Composite layups, material orientations, continuum shells | M8 | yes | planned | -| Model | Element formulation per region (integration, incompatible modes, hybrid) mapped per solver | M3 | yes | planned | -| Model | Boundary conditions: fixed, displacement, symmetry | M3 | | planned | -| Model | Loads: force and moment, pressure, gravity, body force | M3 | | planned | -| Model | Reference points, kinematic and distributing coupling | M3 | | planned | -| Model | Steps: linear static, frequency; output presets | M3 | | planned | +| Model | Element formulation per region (integration, incompatible modes, hybrid) mapped per solver | M3 | yes | done (model; per-solver mapping with the deck writers) | +| Model | Boundary conditions: fixed, displacement, symmetry | M3 | | done | +| Model | Loads: force and moment, pressure, gravity, body force | M3 | | done | +| Model | Reference points, kinematic and distributing coupling | M3 | | done | +| Model | Steps: linear static, frequency; output presets | M3 | | done | | Model | Steps: nonlinear static, explicit dynamics, heat transfer; amplitudes | M7–M8 | | planned | | Interaction | Tie between non-matching meshes | M7 | | planned | | Interaction | Contact (general and pairs), friction, automatic pair detection | M7 | | planned |