diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md
index 83eb313..6912350 100644
--- a/ATTRIBUTION.md
+++ b/ATTRIBUTION.md
@@ -66,6 +66,7 @@ is committed and listed here with its exports.
| assets/icons/mesh.surface.check.svg, assets/icons/mesh.surface.check.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.surface.check.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
| assets/icons/mesh.tag.svg, assets/icons/mesh.tag.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.tag.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
| assets/icons/mesh.tet.generate.svg, assets/icons/mesh.tet.generate.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.tet.generate.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
+| assets/icons/mesh.planar.svg, assets/icons/mesh.planar.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.planar.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
| assets/icons/mesh.transition.svg, assets/icons/mesh.transition.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.transition.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
| assets/icons/mesh.untangle.svg, assets/icons/mesh.untangle.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.untangle.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
| assets/icons/mesh.verify.svg, assets/icons/mesh.verify.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.verify.json` (`cargo xtask icons`) | MIT OR Apache-2.0 |
diff --git a/ROADMAP.md b/ROADMAP.md
index ef21b64..cc908f2 100644
--- a/ROADMAP.md
+++ b/ROADMAP.md
@@ -22,7 +22,7 @@ hex meshes today and it gives the best elements where stresses are read.
| 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: implicit LNC Steiner points on segments, gift-wrapping facet recovery, inside/outside labelling, B-rep face tags on boundary tet faces, few-ulp rounding repair checked by the validity gate; `modelcraft-mesh::cdt` and `mesh.tet.generate`, done: all 20 valid corpus parts); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress |
| M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input (done: `deck.write`), Nastran SOL 101/103 (done: `deck.write`) and LS-DYNA implicit (done); the `.modelcraft` file (text definition, binary mesh, geometry references reopened on load, journal; atomic saves, autosave and recovery; undo/redo on the journal: done) | in progress |
-| M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started |
+| M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion; 2D elements: plane stress, plane strain and axisymmetric tri3, tri6, quad4 (incompatible modes QM6), quad8 with two degrees of freedom per node, centrifugal loads, `mesh.planar` 2D meshes of polygons and planar faces, verified by patch tests, Kirsch, Lamé (plane strain and axisymmetric), rotating discs and cylinders, and a plane stress cantilever; deck mapping to Abaqus CPS/CPE/CAX, Nastran membranes, LS-DYNA ELFORM 12/13); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started |
| 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 |
| M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, springs, point masses; amplitudes, multiple steps; nonlinear static steps in decks | planned |
diff --git a/assets/icons/mesh.planar.svg b/assets/icons/mesh.planar.svg
new file mode 100644
index 0000000..74478ee
--- /dev/null
+++ b/assets/icons/mesh.planar.svg
@@ -0,0 +1,8 @@
+
+
diff --git a/assets/icons/mesh.planar.vectorcraft b/assets/icons/mesh.planar.vectorcraft
new file mode 100644
index 0000000..5bde7ab
Binary files /dev/null and b/assets/icons/mesh.planar.vectorcraft differ
diff --git a/assets/icons/src/mesh.planar.json b/assets/icons/src/mesh.planar.json
new file mode 100644
index 0000000..a22d219
--- /dev/null
+++ b/assets/icons/src/mesh.planar.json
@@ -0,0 +1 @@
+{"shapes": [{"d": "M3 5 L21 5 L21 19 L3 19 Z", "fill": "soft", "stroke": "line"}, {"d": "M9 5 L9 19 M15 5 L15 19 M3 12 L21 12", "stroke": "accent", "width": 1.25}]}
diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs
index 004173b..88a7252 100644
--- a/crates/decks/src/abaqus.rs
+++ b/crates/decks/src/abaqus.rs
@@ -64,6 +64,18 @@
//! | quad8 | | S8R (reduced only) | | |
//! | beam2, rod2 | | | B31 | T3D2 |
//!
+//! | ModelCraft | Plane stress | Plane strain | Axisymmetric |
+//! |---|---|---|---|
+//! | tri3 | CPS3 | CPE3 | CAX3 |
+//! | tri6 | CPS6 | CPE6 | CAX6 |
+//! | quad4 | CPS4, CPS4R, CPS4I | CPE4, CPE4R, CPE4I | CAX4, CAX4R, CAX4I |
+//! | quad8 | CPS8, CPS8R | CPE8, CPE8R | CAX8, CAX8R |
+//!
+//! 2D models: Abaqus's axisymmetric elements, like the model, take r = x and the axis along Y,
+//! and their concentrated loads are totals over the circumference; a 2D section's thickness is
+//! the `*SOLID SECTION` data line. Surfaces of 2D elements are their edges, `S1`… in the edge
+//! order of the model; body forces of axisymmetric elements are `BR`, `BZ`.
+//!
//! `H` marks the hybrid variant. A four-node tetrahedron integrates at one point either way, so
//! reduced integration writes C3D4. Not expressible: reduced integration of ten-node tetrahedra and
//! six-node wedges, full integration of eight-node shells, six-node triangular shells, orthotropic
@@ -125,9 +137,28 @@ pub fn element_type(kind: ElementKind, section: &Section) -> std::result::Result
use SectionKind as S;
let f = section.formulation;
let reduced = f.integration == Integration::Reduced;
- if (f.hybrid || f.incompatible_modes) && section.kind.dimension() != 3 {
+ if (f.hybrid && section.kind.dimension() != 3) || (f.incompatible_modes && !section.kind.is_continuum()) {
return Err(format!("hybrid or incompatible-mode {} sections", section.kind.name()));
}
+ if let Some(prefix) = match section.kind {
+ S::PlaneStress { .. } => Some("CPS"),
+ S::PlaneStrain { .. } => Some("CPE"),
+ S::Axisymmetric {} => Some("CAX"),
+ _ => None,
+ } {
+ let suffix = match (kind, reduced, f.incompatible_modes) {
+ (K::Quad4, false, true) => "4I",
+ (_, true, true) => return Err("incompatible modes with reduced integration".into()),
+ (_, _, true) => return Err(format!("incompatible modes on {} elements (quad4 only)", kind.name())),
+ (K::Tri3, ..) => "3",
+ (K::Tri6, true, _) => return Err("reduced integration of tri6 elements (Abaqus has no reduced-integration variant)".into()),
+ (K::Tri6, ..) => "6",
+ (K::Quad4, r, _) => ["4", "4R"][usize::from(r)],
+ (K::Quad8, r, _) => ["8", "8R"][usize::from(r)],
+ _ => return Err(format!("{} elements in a {} section", kind.name(), section.kind.name())),
+ };
+ return Ok(format!("{prefix}{suffix}"));
+ }
let plain = match (§ion.kind, kind, reduced, f.incompatible_modes) {
(S::Solid {}, K::Hex8, false, true) => "C3D8I",
(S::Solid {}, _, true, true) => return Err("incompatible modes with reduced integration".into()),
@@ -209,7 +240,7 @@ fn check_expressible(model: &Model) -> Result<()> {
}
for (name, section) in model.sections() {
let ortho = matches!(model.material(§ion.material).and_then(|m| m.elastic), Some(Elastic::EngineeringConstants { .. }));
- if ortho && section.kind.dimension() != 3 {
+ if ortho && !section.kind.is_continuum() {
bad.insert(format!(
"section `{name}`: orthotropic material `{}` in a {} section (its material axes need an orientation, not supported yet)",
section.material,
@@ -268,10 +299,22 @@ pub fn face_id(kind: ElementKind, face: u8) -> Option<&'static str> {
}
}
+/// The surface identifier of edge `edge` (1-based) of a 2D element: `S1`… in the model's edge
+/// order (which is Abaqus's face order for plane and axisymmetric elements).
+pub fn edge_id(kind: ElementKind, edge: u8) -> Option<&'static str> {
+ const EDGES: [&str; 4] = ["S1", "S2", "S3", "S4"];
+ let i = usize::from(edge).checked_sub(1)?;
+ if kind.dimension() != 2 || i >= kind.edges().len() {
+ return None;
+ }
+ EDGES.get(i).copied()
+}
+
/// The node surface each coupling ties its reference point to: a fresh name per coupling that
/// clashes with no set (sets with faces own surfaces of their name) and no other coupling surface.
fn coupling_surfaces(model: &Model, r: &Resolved) -> BTreeMap {
- let mut taken: BTreeSet = r.sets.iter().filter(|(_, s)| !s.faces.is_empty()).map(|(n, _)| n.to_ascii_uppercase()).collect();
+ let mut taken: BTreeSet =
+ r.sets.iter().filter(|(n, s)| !s.faces.is_empty() || !r.planar_edges(n).is_empty()).map(|(n, _)| n.to_ascii_uppercase()).collect();
let mut out = BTreeMap::new();
for name in model.couplings().keys() {
let base = format!("{name}-nodes");
@@ -365,6 +408,9 @@ fn element_loads<'a>(loads: &[ElementLoad<'a>], kind: BodyKind, out: &mut Vec<(E
out.retain(|(_, _, v)| v.iter().any(|c| *c != 0.0));
}
+/// A centrifugal load: the element set, ω², a point on the axis and its direction.
+type Centrifugal<'a> = (&'a str, f64, [f64; 3], [f64; 3]);
+
/// The loads of one general step, as the deck writes them.
#[derive(Default)]
struct StepLoads<'a> {
@@ -373,6 +419,9 @@ struct StepLoads<'a> {
/// Surface (set name) → summed pressure.
pressure: BTreeMap<&'a str, f64>,
element: Vec<(ElementTarget<'a>, BodyKind, [f64; 3])>,
+ centrifugal: Vec>,
+ /// Body forces of axisymmetric elements are radial and axial (`BR`, `BZ`).
+ axisymmetric: bool,
}
/// Which load keywords a step wrote entries for, so the next step knows what to remove.
@@ -387,7 +436,8 @@ impl<'a> StepLoads<'a> {
/// `everything` is every element and the section sets that hold them (gravity on the whole
/// model).
fn new(model: &'a Model, r: &'a Resolved, step: &'a Step, everything: &'a (Vec<&'a str>, Vec)) -> StepLoads<'a> {
- let mut s = StepLoads { concentrated: crate::concentrated_loads(model, r, &step.loads), ..StepLoads::default() };
+ let axisymmetric = r.extent() == Some(modelcraft_model::Extent::Revolved);
+ let mut s = StepLoads { concentrated: crate::concentrated_loads(model, r, &step.loads), axisymmetric, ..StepLoads::default() };
let (mut gravity, mut body) = (Vec::new(), Vec::new());
for name in &step.loads {
match model.load(name) {
@@ -397,6 +447,13 @@ impl<'a> StepLoads<'a> {
gravity.push(ElementLoad { elsets: everything.0.clone(), elements: &everything.1, vector: *acceleration })
}
Some(Load::BodyForce { set, vector }) => body.push(ElementLoad::on_set(r, set, *vector)),
+ Some(Load::Centrifugal { angular_velocity, point, axis, set }) => {
+ let w2 = angular_velocity * angular_velocity;
+ match set {
+ Some(set) => s.centrifugal.push((set, w2, *point, *axis)),
+ None => s.centrifugal.extend(everything.0.iter().map(|set| (*set, w2, *point, *axis))),
+ }
+ }
_ => {}
}
}
@@ -407,7 +464,11 @@ impl<'a> StepLoads<'a> {
}
fn written(&self) -> Written {
- Written { concentrated: !self.concentrated.is_empty(), pressure: !self.pressure.is_empty(), element: !self.element.is_empty() }
+ Written {
+ concentrated: !self.concentrated.is_empty(),
+ pressure: !self.pressure.is_empty(),
+ element: !self.element.is_empty() || !self.centrifugal.is_empty(),
+ }
}
}
@@ -513,7 +574,8 @@ impl Deck {
self.line(&format!("*ELSET, ELSET={name}"));
self.labels(&s.elements);
}
- if !s.faces.is_empty() {
+ let edges = r.planar_edges(name);
+ if !s.faces.is_empty() || !edges.is_empty() {
self.line(&format!("*SURFACE, TYPE=ELEMENT, NAME={name}"));
for f in &s.faces {
let id = model.element(f.element).and_then(|e| face_id(e.kind, f.face)).ok_or_else(|| Error::Unsupported {
@@ -522,6 +584,13 @@ impl Deck {
})?;
let _ = writeln!(self.out, "{}, {id}", f.element);
}
+ for g in &edges {
+ let id = model.element(g.element).and_then(|e| edge_id(e.kind, g.edge)).ok_or_else(|| Error::Unsupported {
+ format: FORMAT,
+ what: format!("set `{name}`: element {} has no edge {}", g.element, g.edge),
+ })?;
+ let _ = writeln!(self.out, "{}, {id}", g.element);
+ }
}
}
Ok(())
@@ -565,10 +634,14 @@ impl Deck {
self.comment(&format!("section {name}"));
let (set, material) = (&s.set, &s.material);
match s.kind {
- SectionKind::Solid {} => {
+ SectionKind::Solid {} | SectionKind::Axisymmetric {} => {
self.line(&format!("*SOLID SECTION, ELSET={set}, MATERIAL={material}"));
self.line(",");
}
+ SectionKind::PlaneStress { thickness } | SectionKind::PlaneStrain { thickness } => {
+ self.line(&format!("*SOLID SECTION, ELSET={set}, MATERIAL={material}"));
+ self.numbers_line("", &[thickness])?;
+ }
SectionKind::Rod { area } => {
self.line(&format!("*SOLID SECTION, ELSET={set}, MATERIAL={material}"));
self.numbers_line("", &[area])?;
@@ -686,7 +759,7 @@ impl Deck {
self.numbers_line(&format!("{surface}, P"), &[*p])?;
}
}
- if self.op_new("*DLOAD", !loads.element.is_empty(), before.element) {
+ if self.op_new("*DLOAD", !loads.element.is_empty() || !loads.centrifugal.is_empty(), before.element) {
for (target, kind, v) in &loads.element {
let target = match target {
ElementTarget::Set(s) => s.to_string(),
@@ -698,7 +771,8 @@ impl Deck {
self.numbers_line(&format!("{target}, GRAV"), &[g, v[0] / g, v[1] / g, v[2] / g])?;
}
BodyKind::Force => {
- for (label, c) in ["BX", "BY", "BZ"].into_iter().zip(v) {
+ let labels: &[&str] = if loads.axisymmetric { &["BR", "BZ"] } else { &["BX", "BY", "BZ"] };
+ for (label, c) in labels.iter().zip(v) {
if *c != 0.0 {
self.numbers_line(&format!("{target}, {label}"), &[*c])?;
}
@@ -706,6 +780,10 @@ impl Deck {
}
}
}
+ for (set, w2, point, axis) in &loads.centrifugal {
+ // CENTRIF: ω², a point on the axis and the axis direction.
+ self.numbers_line(&format!("{set}, CENTRIF"), &[*w2, point[0], point[1], point[2], axis[0], axis[1], axis[2]])?;
+ }
}
Ok(())
}
@@ -746,7 +824,7 @@ impl Deck {
self.boundary_lines(model, r, &initial)?;
}
// Gravity on the whole model acts on every element, through the section sets.
- let whole_model_gravity = model.loads().values().any(|l| matches!(l, Load::Gravity { set: None, .. }));
+ let whole_model_gravity = model.loads().values().any(|l| matches!(l, Load::Gravity { set: None, .. } | Load::Centrifugal { set: None, .. }));
let everything: (Vec<&str>, Vec) = if whole_model_gravity {
(model.sections().values().map(|s| s.set.as_str()).collect(), model.elements().keys().copied().collect())
} else {
diff --git a/crates/decks/src/lib.rs b/crates/decks/src/lib.rs
index bbad58e..076d236 100644
--- a/crates/decks/src/lib.rs
+++ b/crates/decks/src/lib.rs
@@ -15,6 +15,7 @@ pub mod nastran;
use std::collections::{BTreeMap, BTreeSet};
+use modelcraft_model::nodal::edge_pressure;
use modelcraft_model::{Diagnostic, Load, Model, Resolved};
#[derive(Debug, Clone, PartialEq, thiserror::Error)]
@@ -134,3 +135,26 @@ pub fn concentrated_loads<'a>(model: &Model, r: &Resolved, active: impl IntoIter
out.retain(|_, v| *v != 0.0);
out
}
+
+/// The consistent nodal forces (degrees of freedom 1 and 2) of the active pressures on a 2D
+/// model's element edges, for formats that write them as concentrated loads: the model's edge
+/// integration ([`modelcraft_model::nodal::edge_pressure`]), a positive pressure pushing
+/// into the element, summed per node and direction; zero sums are left out.
+pub fn edge_pressure_loads<'a>(model: &Model, r: &Resolved, active: impl IntoIterator- ) -> Result> {
+ let mut out: BTreeMap<(u32, u8), f64> = BTreeMap::new();
+ for name in active {
+ let Some(Load::Pressure { set, magnitude }) = model.load(name) else { continue };
+ for g in r.planar_edges(set) {
+ let Some(extent) = r.planar.get(&g.element).copied() else { continue };
+ let labels = model.edge_nodes(g).unwrap_or_default();
+ let xyz: Vec<[f64; 3]> = labels.iter().filter_map(|n| model.node(*n)).map(|n| n.xyz).collect();
+ for (node, f) in labels.iter().zip(edge_pressure(&xyz, extent, *magnitude)?) {
+ for (dof, c) in (1u8..).zip(&f[..2]) {
+ *out.entry((*node, dof)).or_default() += c;
+ }
+ }
+ }
+ }
+ out.retain(|_, v| *v != 0.0);
+ Ok(out)
+}
diff --git a/crates/decks/src/lsdyna/elements.rs b/crates/decks/src/lsdyna/elements.rs
index 11b10bd..980dbc6 100644
--- a/crates/decks/src/lsdyna/elements.rs
+++ b/crates/decks/src/lsdyna/elements.rs
@@ -97,7 +97,9 @@ pub(crate) struct Elform {
/// hexahedra), reduced 15 (two-point pentahedron); tet4 10 (one point, exact for a linear
/// tetrahedron), hybrid 13 (nodal pressure); tet10 16; hex20 23. Shells: quad4 16 (fully
/// integrated), reduced 2 (one-point) with hourglass control; tri3 4 (C0 triangle); quad8 23;
-/// tri6 24. Rods: 3 (truss).
+/// tri6 24. Rods: 3 (truss). 2D (in the XY plane): plane stress 12, plane strain 13, both one-point
+/// forms on quad4 (written for reduced integration, with hourglass control) and tri3 (exact at one
+/// point); quadratic 2D elements and axisymmetric sections have no form written here.
/// What it cannot express is the reason, naming the section.
pub(crate) fn elform(kind: ElementKind, section: &SectionKind, f: Formulation, section_name: &str) -> Result {
let plain = |elform| Ok(Elform { elform, hourglass: None });
@@ -142,6 +144,18 @@ pub(crate) fn elform(kind: ElementKind, section: &SectionKind, f: Formulation, s
_ => no("a rod section cannot hold"),
},
SectionKind::Beam { .. } => no("beam sections are not written to LS-DYNA decks yet (planned with beam profiles):"),
+ SectionKind::Axisymmetric {} => {
+ no("axisymmetric sections are not written to LS-DYNA decks yet (its axisymmetric forms count loads per radian; not verified):")
+ }
+ SectionKind::PlaneStress { .. } | SectionKind::PlaneStrain { .. } => {
+ let form = if matches!(section, SectionKind::PlaneStress { .. }) { 12 } else { 13 };
+ match kind {
+ ElementKind::Tri3 => plain(form),
+ ElementKind::Quad4 if reduced => Ok(Elform { elform: form, hourglass: Some(Hourglass::Shell) }),
+ ElementKind::Quad4 => no("LS-DYNA's 2D forms integrate at one point (set the section's integration to reduced) for"),
+ _ => no("LS-DYNA's 2D forms are linear; there is no form for"),
+ }
+ }
}
}
diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs
index 5c27750..8721b34 100644
--- a/crates/decks/src/lsdyna/mod.rs
+++ b/crates/decks/src/lsdyna/mod.rs
@@ -30,6 +30,11 @@
//! order when LS-DYNA is started with `CASE` on its command line. A single step is written
//! without cases.
//!
+//! 2D models (plane stress and plane strain in the XY plane) are `*SECTION_SHELL` parts with
+//! ELFORM 12 and 13 ([`elements`]); plane strain is per unit thickness in LS-DYNA, so only a
+//! thickness of 1 is written; pressures on 2D element edges are the consistent nodal forces
+//! (`*LOAD_NODE_POINT`).
+//!
//! What LS-DYNA cannot take is an [`Error::Unsupported`] naming the item: beam sections (not yet
//! written), body forces per volume, kinematic couplings of some degrees of freedom only, boundary
//! conditions on a coupled reference point, gravity on different sets, hybrid formulations other
@@ -188,6 +193,11 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result> {
let el = model.element(*e).ok_or_else(|| broken(format!("element {e} of section `{name}`")))?;
by_kind.entry(el.kind).or_default().push((*e, el));
}
+ if let SectionKind::PlaneStrain { thickness } = sec.kind
+ && thickness != 1.0
+ {
+ refused.insert(format!("section `{name}`: a plane strain thickness of {thickness} (LS-DYNA's plane strain is per unit thickness)"));
+ }
if let SectionKind::Shell { points, .. } = sec.kind
&& points > MAX_SHELL_POINTS
{
@@ -270,9 +280,13 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result> {
}
}
Load::Traction { .. } => {}
+ Load::Pressure { .. } if r.is_planar() => {}
Load::Pressure { set, .. } => {
segment_sets.insert(set.as_str());
}
+ Load::Centrifugal { .. } => {
+ refused.insert(format!("load `{l}`: a centrifugal load (not written to LS-DYNA decks yet)"));
+ }
Load::Gravity { set, .. } => {
gravity_sets.insert(set.as_deref());
}
@@ -414,6 +428,14 @@ fn parts(d: &mut Deck, plan: &Plan) -> Result<()> {
d.card([Field::from(area)]);
}
SectionKind::Beam { .. } => return Err(broken(format!("section `{}`: beams have no LS-DYNA formulation", p.name))),
+ SectionKind::PlaneStress { thickness } | SectionKind::PlaneStrain { thickness } => {
+ d.line("*SECTION_SHELL");
+ d.header(&["secid", "elform"]);
+ d.card([Field::from(p.id), Field::from(p.elform)]);
+ d.header(&["t1", "t2", "t3", "t4"]);
+ d.card([Field::from(thickness), Field::from(thickness), Field::from(thickness), Field::from(thickness)]);
+ }
+ SectionKind::Axisymmetric {} => return Err(broken(format!("section `{}`: axisymmetric sections have no LS-DYNA formulation", p.name))),
}
}
for (h, hgid) in &plan.hourglasses {
@@ -621,12 +643,14 @@ fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Ste
if is_static {
let mut gravity = [0.0f64; 3];
let mut pointwise = Vec::new();
+ let mut edge_pressures = Vec::new();
for (name, load) in &loads {
match load {
Load::Force { .. } | Load::Moment { .. } if plan.uniform.contains_key(name.as_str()) => {
on_node_set(d, plan, name, matches!(load, Load::Moment { .. }), n);
}
Load::Force { .. } | Load::Moment { .. } | Load::Traction { .. } => pointwise.push(*name),
+ Load::Pressure { .. } if r.is_planar() => edge_pressures.push(name.as_str()),
Load::Pressure { set, magnitude } => {
let sid = plan.sets.get(set.as_str()).copied().unwrap_or_default();
d.comment(&format!("load `{name}`"));
@@ -640,10 +664,15 @@ fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Ste
}
}
Load::BodyForce { .. } => return Err(broken(format!("load `{name}`: a body force"))),
+ Load::Centrifugal { .. } => return Err(broken(format!("load `{name}`: a centrifugal load"))),
}
}
// The other concentrated loads node by node, summed where they share a node.
- let point = crate::concentrated_loads(model, r, pointwise);
+ let mut point = crate::concentrated_loads(model, r, pointwise);
+ for (key, v) in crate::edge_pressure_loads(model, r, edge_pressures)? {
+ *point.entry(key).or_default() += v;
+ }
+ point.retain(|_, v| *v != 0.0);
if !point.is_empty() {
d.line("*LOAD_NODE_POINT");
d.header(&["nid", "dof", "lcid", "sf"]);
diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs
index 5b86fb9..be07584 100644
--- a/crates/decks/src/nastran/mod.rs
+++ b/crates/decks/src/nastran/mod.rs
@@ -24,7 +24,15 @@
//! gravity: `GRAV`; a step with several loads combines them with `LOAD`;
//! - kinematic coupling: `RBE2` (the reference point independent); distributing coupling: `RBE3`
//! (equal weights on the set's translations; the reference point dependent);
-//! - frequency steps: `EIGRL` (Lanczos) with the step's mode count and range.
+//! - frequency steps: `EIGRL` (Lanczos) with the step's mode count and range;
+//! - 2D models: plane stress as membranes, `CTRIA3`, `CQUAD4`, `CTRIA6`, `CQUAD8` with a
+//! `PSHELL` that has only a membrane material, every grid holding its out-of-plane degrees of
+//! freedom 3–6 (`GRID` field PS); edge pressures as the consistent `FORCE`s the model's edge
+//! integration gives. Plane strain and axisymmetric elements are not written: their cards
+//! (`PPLANE`/`PLPLANE` with `CPLSTN*`/`CQUADX`, `CTRIAX6`, `PAXI`) differ between Nastran
+//! dialects, so a deck could not be read the same everywhere;
+//! - centrifugal loads: `RFORCE` about an axis through the origin of the basic system (speed in
+//! revolutions per time unit).
//!
//! Field widths, continuations and number forms are in [`card`]. Model content Nastran cannot
//! express is reported as [`Error::Unsupported`], every item at once, and nothing is written.
@@ -188,10 +196,31 @@ struct MaterialUse {
/// One property entry, and what its elements' connection entries need besides their nodes.
enum Property {
- Solid { mid: u32, rule: (Option<&'static str>, Option<&'static str>) },
- Shell { mid: u32, thickness: f64, offset: Option, mcid: Option },
- Beam { mid: u32, profile: BeamProfile, orientation: [f64; 3] },
- Rod { mid: u32, area: f64 },
+ Solid {
+ mid: u32,
+ rule: (Option<&'static str>, Option<&'static str>),
+ },
+ Shell {
+ mid: u32,
+ thickness: f64,
+ offset: Option,
+ mcid: Option,
+ },
+ Beam {
+ mid: u32,
+ profile: BeamProfile,
+ orientation: [f64; 3],
+ },
+ Rod {
+ mid: u32,
+ area: f64,
+ },
+ /// Plane stress: a shell with only a membrane material.
+ Membrane {
+ mid: u32,
+ thickness: f64,
+ mcid: Option,
+ },
}
/// A constraint set: per node, the prescribed value of each degree of freedom.
@@ -299,7 +328,9 @@ impl<'m> Deck<'m> {
let u = uses.entry(s.material.as_str()).or_default();
match s.kind {
SectionKind::Solid {} => u.solid = true,
- SectionKind::Shell { .. } => u.shell = true,
+ SectionKind::Shell { .. } | SectionKind::PlaneStress { .. } | SectionKind::PlaneStrain { .. } | SectionKind::Axisymmetric {} => {
+ u.shell = true
+ }
SectionKind::Beam { .. } | SectionKind::Rod { .. } => u.line = true,
}
}
@@ -398,6 +429,25 @@ impl<'m> Deck<'m> {
let Some(mid) = ids.isotropic else { continue };
self.add_property(sname, Property::Rod { mid, area }, elements.iter().copied())?;
}
+ SectionKind::PlaneStress { thickness } => {
+ if sec.formulation.integration == Integration::Reduced || sec.formulation.incompatible_modes {
+ self.problems.insert(format!(
+ "section `{sname}`: reduced integration or incompatible modes of plane stress elements (Nastran's membranes have one formulation)"
+ ));
+ }
+ let (mid, mcid) = match ids {
+ MaterialIds { isotropic: Some(m), .. } => (m, None),
+ MaterialIds { shell: Some(m), .. } => (m, Some(*self.material_csys.get_or_insert(1))),
+ _ => continue,
+ };
+ self.add_property(sname, Property::Membrane { mid, thickness, mcid }, elements.iter().copied())?;
+ }
+ SectionKind::PlaneStrain { .. } | SectionKind::Axisymmetric {} => {
+ self.problems.insert(format!(
+ "section `{sname}`: {} elements (Nastran's plane strain and axisymmetric cards differ between dialects; not written)",
+ sec.kind.name()
+ ));
+ }
}
}
Ok(())
@@ -529,6 +579,22 @@ impl<'m> Deck<'m> {
self.problems.insert(format!("load `{name}`: gravity on part of the model (set `{s}`); Nastran's gravity acts on all of it"));
}
Load::Gravity { .. } => gravity += 1,
+ Load::Centrifugal { set: Some(s), .. } if r.sets.get(s).map(|s| s.elements.len()) != Some(model.elements().len()) => {
+ self.problems
+ .insert(format!("load `{name}`: a rotation of part of the model (set `{s}`); Nastran's RFORCE acts on all of it"));
+ }
+ Load::Centrifugal { point, axis, .. } => {
+ let c = [
+ point[1] * axis[2] - point[2] * axis[1],
+ point[2] * axis[0] - point[0] * axis[2],
+ point[0] * axis[1] - point[1] * axis[0],
+ ];
+ if c.iter().any(|v| *v != 0.0) {
+ self.problems.insert(format!(
+ "load `{name}`: a rotation about an axis that misses the origin (written as RFORCE about the origin)"
+ ));
+ }
+ }
_ => {}
}
let sid = next_id(&mut self.next_sid)?;
@@ -597,7 +663,10 @@ impl<'m> Deck<'m> {
fn emit_mesh(&self, out: &mut String) -> Result<()> {
comment(out, "nodes");
for (label, node) in self.model.nodes() {
- Card::new("GRID").int(*label).blank().reals(node.xyz).write(out)?;
+ let grid = Card::new("GRID").int(*label).blank().reals(node.xyz);
+ // A 2D model's grids hold their out-of-plane translation and rotations.
+ let grid = if self.r.is_planar() { grid.blank().int(3456) } else { grid };
+ grid.write(out)?;
}
if !self.model.reference_points().is_empty() {
comment(out, "reference points");
@@ -650,6 +719,7 @@ impl<'m> Deck<'m> {
Card::new("PBEAML").int(pid).int(*mid).blank().text(shape).end_line().reals(dims)
}
Property::Rod { mid, area } => Card::new("PROD").int(pid).int(*mid).real(*area),
+ Property::Membrane { mid, thickness, .. } => Card::new("PSHELL").int(pid).int(*mid).real(*thickness),
};
write_named(out, &format!("section {section}"), pid, card)?;
}
@@ -664,6 +734,8 @@ impl<'m> Deck<'m> {
// CQUAD8 has four corner thicknesses before the material axes and the offset.
Property::Shell { offset, mcid, .. } if el.kind == ElementKind::Quad8 => card.blank().blank().blank().blank().opt(*mcid).opt(*offset),
Property::Shell { offset, mcid, .. } => card.opt(*mcid).opt(*offset),
+ Property::Membrane { mcid, .. } if el.kind == ElementKind::Quad8 => card.blank().blank().blank().blank().opt(*mcid),
+ Property::Membrane { mcid, .. } => card.opt(*mcid),
Property::Solid { .. } | Property::Rod { .. } => card,
};
card.write(out)?;
@@ -748,6 +820,18 @@ impl<'m> Deck<'m> {
entry(Card::new(card).int(*sid).int(*node).blank().real(1.0).reals(*v), out, &mut first)?;
}
}
+ Load::Pressure { .. } if self.r.is_planar() => {
+ let forces = crate::edge_pressure_loads(self.model, self.r, [*name])?;
+ let mut by_node: BTreeMap = BTreeMap::new();
+ for ((node, dof), v) in forces {
+ if let Some(c) = by_node.entry(node).or_default().get_mut(usize::from(dof) - 1) {
+ *c = v;
+ }
+ }
+ for (node, v) in by_node {
+ entry(Card::new("FORCE").int(*sid).int(node).blank().real(1.0).reals(v), out, &mut first)?;
+ }
+ }
Load::Pressure { set, magnitude } => {
for f in self.r.sets.get(set).map(|s| s.faces.as_slice()).unwrap_or_default() {
if let Some(card) = self.pressure(*sid, *f, *magnitude) {
@@ -756,6 +840,12 @@ impl<'m> Deck<'m> {
}
}
Load::Gravity { acceleration, .. } => entry(Card::new("GRAV").int(*sid).blank().real(1.0).reals(*acceleration), out, &mut first)?,
+ Load::Centrifugal { angular_velocity, axis, .. } => {
+ // RFORCE SID G CID A R1 R2 R3: about the basic origin (G blank), A in
+ // revolutions per time unit along the direction R.
+ let a = angular_velocity / (2.0 * std::f64::consts::PI);
+ entry(Card::new("RFORCE").int(*sid).blank().blank().real(a).reals(*axis), out, &mut first)?
+ }
Load::BodyForce { .. } => {}
}
}
diff --git a/crates/decks/tests/planar.rs b/crates/decks/tests/planar.rs
new file mode 100644
index 0000000..01bd09f
--- /dev/null
+++ b/crates/decks/tests/planar.rs
@@ -0,0 +1,157 @@
+//! 2D models (plane stress, plane strain, axisymmetric) in every deck format: element names,
+//! section data, edge surfaces and pressures, centrifugal loads, and what a format cannot write.
+
+use modelcraft_decks::{Error, abaqus, edge_pressure_loads, lsdyna, nastran};
+use modelcraft_model::{
+ BoundaryCondition, Elastic, ElementEdge, ElementKind, Formulation, Integration, Load, Material, Model, Procedure, Section, SectionKind, Set,
+ Step, Target, UnitSystem,
+};
+
+/// Two quad4 over [1, 3] × [0, 1] (or the given kind on the same corners), the left edge held,
+/// a pressure on the right edge and a centrifugal load.
+fn plate(kind: SectionKind, element: ElementKind, formulation: Formulation) -> Model {
+ let mut m = Model::default();
+ m.set_units(UnitSystem::MmTonneS);
+ for (i, x) in [1.0, 2.0, 3.0].iter().enumerate() {
+ m.set_node(i as u32 + 1, [*x, 0.0, 0.0]).unwrap();
+ m.set_node(i as u32 + 4, [*x, 1.0, 0.0]).unwrap();
+ }
+ match element {
+ ElementKind::Tri3 => {
+ for (l, n) in [(1, [1, 2, 5]), (2, [1, 5, 4]), (3, [2, 3, 6]), (4, [2, 6, 5])] {
+ m.set_element(l, ElementKind::Tri3, n.to_vec()).unwrap();
+ }
+ }
+ _ => {
+ m.set_element(1, ElementKind::Quad4, vec![1, 2, 5, 4]).unwrap();
+ m.set_element(2, ElementKind::Quad4, vec![2, 3, 6, 5]).unwrap();
+ }
+ }
+ let right = if element == ElementKind::Tri3 { ElementEdge { element: 3, edge: 2 } } else { ElementEdge { element: 2, edge: 2 } };
+ let all: Vec = m.elements().keys().copied().collect();
+ m.add_set("all", Set { elements: all, ..Set::default() }).unwrap();
+ m.add_set("left", Set { nodes: vec![1, 4], ..Set::default() }).unwrap();
+ m.add_set("right", Set { edges: vec![right], ..Set::default() }).unwrap();
+ m.add_material("Steel", Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), fluid: None })
+ .unwrap();
+ m.add_section("plate", Section { set: "all".into(), material: "Steel".into(), kind, formulation }).unwrap();
+ m.add_boundary_condition("hold", BoundaryCondition::Fixed { target: Target::Set("left".into()) }).unwrap();
+ m.add_load("press", Load::Pressure { set: "right".into(), magnitude: 5.0 }).unwrap();
+ let axis = if kind == (SectionKind::Axisymmetric {}) { [0.0, 1.0, 0.0] } else { [0.0, 0.0, 1.0] };
+ m.add_load("spin", Load::Centrifugal { angular_velocity: 100.0, point: [0.0; 3], axis, set: None }).unwrap();
+ let mut s = Step::new("Load", Procedure::Static {});
+ s.boundary_conditions.insert("hold".into());
+ s.loads.extend(["press".to_string(), "spin".to_string()]);
+ m.add_step(s, None).unwrap();
+ m
+}
+
+fn unsupported(r: Result) -> String {
+ match r {
+ Err(Error::Unsupported { what, .. }) => what,
+ other => panic!("expected unsupported, got {other:?}"),
+ }
+}
+
+#[test]
+fn abaqus_writes_plane_and_axisymmetric_elements() {
+ let full = Formulation::default();
+ let im = Formulation { incompatible_modes: true, ..Formulation::default() };
+ let red = Formulation { integration: Integration::Reduced, ..Formulation::default() };
+ for (kind, f, want) in [
+ (SectionKind::PlaneStress { thickness: 2.0 }, full, "CPS4"),
+ (SectionKind::PlaneStress { thickness: 2.0 }, im, "CPS4I"),
+ (SectionKind::PlaneStrain { thickness: 1.0 }, red, "CPE4R"),
+ (SectionKind::Axisymmetric {}, im, "CAX4I"),
+ ] {
+ let deck = abaqus::write(&plate(kind, ElementKind::Quad4, f)).unwrap();
+ assert!(deck.contains(&format!("*ELEMENT, TYPE={want}\n")), "{want}\n{deck}");
+ // The right edge is surface S2 of element 2; the pressure is a surface load on it.
+ assert!(deck.contains("*SURFACE, TYPE=ELEMENT, NAME=right\n2, S2\n"), "{deck}");
+ assert!(deck.contains("right, P, 5.0\n"), "{deck}");
+ assert!(deck.contains("all, CENTRIF, 10000.0, 0.0, 0.0, 0.0, "), "{deck}");
+ let section = match kind {
+ SectionKind::Axisymmetric {} => "*SOLID SECTION, ELSET=all, MATERIAL=Steel\n,\n",
+ SectionKind::PlaneStress { .. } => "*SOLID SECTION, ELSET=all, MATERIAL=Steel\n2.0\n",
+ _ => "*SOLID SECTION, ELSET=all, MATERIAL=Steel\n1.0\n",
+ };
+ assert!(deck.contains(section), "{deck}");
+ }
+ let deck = abaqus::write(&plate(SectionKind::PlaneStrain { thickness: 1.0 }, ElementKind::Tri3, full)).unwrap();
+ assert!(deck.contains("*ELEMENT, TYPE=CPE3\n") && deck.contains("3, S2\n"), "{deck}");
+ // Axisymmetric body forces are radial and axial.
+ let mut m = plate(SectionKind::Axisymmetric {}, ElementKind::Quad4, full);
+ m.add_load("swell", Load::BodyForce { set: "all".into(), vector: [1.0, 2.0, 0.0] }).unwrap();
+ m.set_load_active("Load", "swell", true).unwrap();
+ let deck = abaqus::write(&m).unwrap();
+ assert!(deck.contains("all, BR, 1.0\nall, BZ, 2.0\n"), "{deck}");
+ // Abaqus has no reduced six-node triangle.
+ let mut m = plate(SectionKind::PlaneStress { thickness: 1.0 }, ElementKind::Quad4, red);
+ for (l, n) in [(10, 10), (11, 11), (12, 12), (13, 13), (14, 14), (15, 15)] {
+ m.set_node(l, [1.0 + f64::from(n - 10) * 0.1, 2.0 + f64::from(n % 2), 0.0]).unwrap();
+ }
+ m.set_element(9, ElementKind::Tri6, vec![10, 12, 14, 11, 13, 15]).unwrap();
+ m.replace_set("all", Set { elements: vec![1, 2, 9], ..Set::default() }).unwrap();
+ assert!(abaqus::write(&m).is_err());
+}
+
+#[test]
+fn nastran_writes_plane_stress_membranes_and_refuses_the_rest() {
+ let m = plate(SectionKind::PlaneStress { thickness: 2.0 }, ElementKind::Quad4, Formulation::default());
+ let deck = nastran::write(&m).unwrap();
+ // A membrane: the material on MID1 only.
+ assert!(deck.contains("PSHELL* 1 1 2.\n$"), "{deck}");
+ assert!(deck.contains("CQUAD4 1 1 1 2 5 4"), "{deck}");
+ // Every grid holds its out-of-plane degrees of freedom (PS on the continuation).
+ let grids = deck.lines().filter(|l| l.starts_with("GRID")).count();
+ assert_eq!(grids, 6);
+ assert_eq!(deck.lines().filter(|l| l.starts_with('*') && l.trim_end().ends_with("3456")).count(), grids, "{deck}");
+ // The edge pressure as consistent forces, pushing in (−X on the right edge).
+ let forces: Vec<&str> = deck.lines().filter(|l| l.starts_with("FORCE")).collect();
+ assert_eq!(forces.len(), 2, "{deck}");
+ assert!(deck.contains("RFORCE"), "{deck}");
+ for kind in [SectionKind::PlaneStrain { thickness: 1.0 }, SectionKind::Axisymmetric {}] {
+ let why = unsupported(nastran::write(&plate(kind, ElementKind::Quad4, Formulation::default())));
+ assert!(why.contains("differ between dialects"), "{why}");
+ }
+ let mut m = plate(SectionKind::PlaneStress { thickness: 2.0 }, ElementKind::Quad4, Formulation::default());
+ m.replace_load("spin", Load::Centrifugal { angular_velocity: 1.0, point: [5.0, 0.0, 0.0], axis: [0.0, 0.0, 1.0], set: None }).unwrap();
+ assert!(unsupported(nastran::write(&m)).contains("misses the origin"));
+}
+
+#[test]
+fn lsdyna_writes_one_point_plane_forms_and_refuses_the_rest() {
+ let red = Formulation { integration: Integration::Reduced, ..Formulation::default() };
+ let mut m = plate(SectionKind::PlaneStress { thickness: 2.0 }, ElementKind::Quad4, red);
+ m.remove_step("Load").unwrap();
+ let mut s = Step::new("Load", Procedure::Static {});
+ s.boundary_conditions.insert("hold".into());
+ s.loads.insert("press".into());
+ m.add_step(s, None).unwrap();
+ let deck = lsdyna::write(&m).unwrap().text;
+ assert!(deck.contains("*SECTION_SHELL") && deck.contains("*HOURGLASS"), "{deck}");
+ assert!(deck.contains("*LOAD_NODE_POINT") && !deck.contains("*LOAD_SEGMENT"), "{deck}");
+ let why = |m: &Model| match lsdyna::write(m) {
+ Err(Error::Unsupported { what, .. }) => what,
+ other => panic!("{other:?}"),
+ };
+ assert!(why(&plate(SectionKind::PlaneStress { thickness: 2.0 }, ElementKind::Quad4, Formulation::default())).contains("one point"));
+ assert!(why(&plate(SectionKind::Axisymmetric {}, ElementKind::Tri3, red)).contains("axisymmetric"));
+ assert!(why(&plate(SectionKind::PlaneStrain { thickness: 2.0 }, ElementKind::Tri3, red)).contains("per unit thickness"));
+ assert!(why(&plate(SectionKind::PlaneStrain { thickness: 1.0 }, ElementKind::Tri3, red)).contains("centrifugal"));
+}
+
+#[test]
+fn edge_pressures_are_the_consistent_nodal_forces() {
+ // Pressure 5 on the right edge (x = 3, length 1): plane stress t = 2 gives −5 per node in X;
+ // axisymmetric, the ring 2π · 3 · 1 shares 5 · 6π equally on a straight edge at fixed r.
+ let m = plate(SectionKind::PlaneStress { thickness: 2.0 }, ElementKind::Quad4, Formulation::default());
+ let r = m.resolve().unwrap();
+ let f = edge_pressure_loads(&m, &r, ["press"]).unwrap();
+ assert_eq!(f.len(), 2);
+ assert!(f.values().all(|v| (v + 5.0).abs() < 1e-12), "{f:?}");
+ let m = plate(SectionKind::Axisymmetric {}, ElementKind::Quad4, Formulation::default());
+ let r = m.resolve().unwrap();
+ let total: f64 = edge_pressure_loads(&m, &r, ["press"]).unwrap().values().sum();
+ assert!((total + 5.0 * 6.0 * std::f64::consts::PI).abs() < 1e-10, "{total}");
+}
diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs
index c9fb222..b71ede9 100644
--- a/crates/engine/src/catalogue.rs
+++ b/crates/engine/src/catalogue.rs
@@ -168,6 +168,13 @@ pub static CATALOGUE: &[Capability] = &[
),
soon("mesh.order.quadratic", "Quadratic", Mesh, "Tet", "M2", "Second-order elements with mid-side nodes on the geometry"),
ready("mesh.delaunay.points", "Delaunay Points", Mesh, "Tet", "Delaunay tetrahedralization of a point list (debugging aid)"),
+ ready(
+ "mesh.planar",
+ "2D Mesh",
+ Mesh,
+ "2D",
+ "Quad or triangle mesh of a planar face or polygon in the XY plane for plane stress, plane strain and axisymmetric models; edges tagged",
+ ),
soon("mesh.quad.faces", "Quad Faces", Mesh, "Structured Hex", "M5", "All-quad meshes of faces on shared edge nodes"),
soon("mesh.hex.mapped", "Mapped Hex", Mesh, "Structured Hex", "M5", "Mapped hex mesh of six-sided regions"),
soon("mesh.hex.sweep", "Sweep Hex", Mesh, "Structured Hex", "M5", "Hex mesh swept along extrusions and revolutions"),
@@ -213,7 +220,7 @@ pub static CATALOGUE: &[Capability] = &[
"Section",
Properties,
"Sections",
- "Solid, shell, beam or rod section assigned to a set, with its element formulation",
+ "Solid, shell, beam, rod or 2D (plane stress, plane strain, axisymmetric) section assigned to a set, with its element formulation",
&["model.section.edit", "model.section.delete"],
),
soon("model.section.composite", "Composite", Properties, "Sections", "M8", "Composite layups and continuum shells"),
diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs
index 8541962..4c6d09c 100644
--- a/crates/engine/src/cmd/define.rs
+++ b/crates/engine/src/cmd/define.rs
@@ -66,7 +66,7 @@ definition_commands!(
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`, …)"#,
+ 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}], "edges": [{"element": 1, "edge": 2}]} — any of the member lists (element edges: where a 2D model's surface loads act); 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,
@@ -94,7 +94,7 @@ definition_commands!(
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}}"#,
+ 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} | {"type": "plane_stress", "thickness": 1} | {"type": "plane_strain", "thickness": 1 (default)} | {"type": "axisymmetric"} (2D: elements in the XY plane; axisymmetric about Y with x the radius, forces totals over the circumference), "formulation": {"integration": "full|reduced", "incompatible_modes": false, "hybrid": false}}"#,
add_section,
replace_section,
remove_section,
@@ -150,7 +150,7 @@ definition_commands!(
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)"#,
+ 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)} | {"name", "type": "centrifugal", "angular_velocity" (rad per time unit), "point": [x, y, z] (on the axis), "axis": [x, y, z], "set" (optional)} — 2D models: in-plane vectors, surface loads on element edges; active in the steps that list it (model.step.activate)"#,
add_load,
replace_load,
remove_load,
diff --git a/crates/engine/src/cmd/mesh.rs b/crates/engine/src/cmd/mesh.rs
index 8235a60..2bd8ef0 100644
--- a/crates/engine/src/cmd/mesh.rs
+++ b/crates/engine/src/cmd/mesh.rs
@@ -1,4 +1,5 @@
use std::borrow::Cow;
+use std::collections::BTreeMap;
use modelcraft_geometry::SurfaceMesh;
use modelcraft_mesh::cdt::{CdtOptions, tetrahedralize};
@@ -87,6 +88,17 @@ pub const TET_GENERATE: CommandSpec = CommandSpec {
run: tet_generate,
};
+/// The first node and element labels of a mesh of `nodes` nodes and `elements` elements added
+/// after the model's largest labels, or an error when the labels would pass [`MAX_LABEL`].
+fn next_labels(s: &Session, nodes: usize, elements: usize) -> Result<(u32, u32)> {
+ let next = |last: Option| last.map_or(Some(1), |l| l.checked_add(1));
+ let fits = |first: Option, count: usize| first.filter(|&f| u64::from(f) + count as u64 - 1 <= u64::from(MAX_LABEL));
+ match (fits(next(s.model.nodes().keys().next_back().copied()), nodes), fits(next(s.model.elements().keys().next_back().copied()), elements)) {
+ (Some(n), Some(e)) => Ok((n, e)),
+ _ => Err(EngineError::Params(format!("the mesh needs labels beyond {MAX_LABEL}"))),
+ }
+}
+
/// The most Steiner points and tetrahedra `mesh.tet.generate` may be asked for.
const TET_MAX_STEINER: u64 = 10_000_000;
const TET_MAX_TETS: u64 = 50_000_000;
@@ -127,13 +139,7 @@ fn tet_generate(s: &mut Session, p: &Value) -> Result {
let cdt = tetrahedralize(&surface, &cdt_options)?;
let gate = cdt.gate(surface.volume())?.into_result()?;
- let next = |last: Option| last.map_or(Some(1), |l| l.checked_add(1));
- let first_node = next(s.model.nodes().keys().next_back().copied());
- let first_element = next(s.model.elements().keys().next_back().copied());
- let fits = |first: Option, count: usize| first.filter(|&f| u64::from(f) + count as u64 - 1 <= u64::from(MAX_LABEL));
- let (Some(first_node), Some(first_element)) = (fits(first_node, cdt.points.len()), fits(first_element, cdt.tets.len())) else {
- return Err(EngineError::Params(format!("the mesh needs labels beyond {MAX_LABEL}")));
- };
+ let (first_node, first_element) = next_labels(s, cdt.points.len(), cdt.tets.len())?;
for (i, xyz) in cdt.points.iter().enumerate() {
s.model.set_node(first_node + i as u32, *xyz)?;
}
@@ -175,6 +181,143 @@ fn tet_generate(s: &mut Session, p: &Value) -> Result {
}))
}
+pub const PLANAR: CommandSpec = CommandSpec {
+ id: "mesh.planar",
+ label: "2D Mesh",
+ params: r#"{"size": 1.0, "kind": "quad4|quad8|tri3|tri6" (default quad4), and either "polygon": [[x, y], ...] with optional "holes": [[[x, y], ...], ...] and "name": "Plane" (the body name its tags use), or "body": "Body1", "face": "F2:top", "tolerance": 0.05 (a planar face of an open design, parallel to XY); optional "max_elements": 1000000 (1..=4000000)} — meshes the planar region with 2D elements in the XY plane (a face at z = c is meshed at z = 0) for plane stress, plane strain or axisymmetric sections (x the radius, Y the axis): free quadrilaterals (or constrained Delaunay triangles) of about `size`, mid-side nodes on the boundary for quadratic kinds. Adds nodes and elements after the largest existing labels and tags them: the elements with the face (polygon: face `area` of body `name`), the boundary element edges with edges named after the neighbouring face of the body (polygon: `e1`, `e2`, … side by side, the outer loop then each hole, and vertices `v1`, … at the corners), so sets on those edges carry loads and boundary conditions"#,
+ kind: CommandKind::Model,
+ run: planar,
+};
+
+/// Most points of a `mesh.planar` polygon or hole, and most holes.
+const PLANAR_MAX_POINTS: usize = 100_000;
+const PLANAR_MAX_HOLES: usize = 10_000;
+const PLANAR_MAX_ELEMENTS: u64 = 4_000_000;
+
+/// A list of `[x, y]` points.
+fn polygon(v: &Value, what: &str) -> Result> {
+ let bad = || EngineError::Params(format!("`{what}` must be a list of 3 to {PLANAR_MAX_POINTS} [x, y] points"));
+ let raw = v.as_array().filter(|a| (3..=PLANAR_MAX_POINTS).contains(&a.len())).ok_or_else(bad)?;
+ raw.iter()
+ .map(|p| match p.as_array().map(Vec::as_slice) {
+ Some([x, y]) => x.as_f64().zip(y.as_f64()).map(|(x, y)| [x, y]).filter(|c| c.iter().all(|v| v.is_finite())),
+ _ => None,
+ })
+ .collect::>>()
+ .ok_or_else(bad)
+}
+
+fn planar(s: &mut Session, p: &Value) -> Result {
+ use modelcraft_mesh::planar::{Curve, face_boundary, mesh_region};
+ use modelcraft_model::Entity;
+
+ let size = p
+ .get("size")
+ .and_then(Value::as_f64)
+ .filter(|v| *v > 0.0 && v.is_finite())
+ .ok_or_else(|| EngineError::Params("`size` must be a positive number (the element size)".into()))?;
+ let kind = match p.get("kind") {
+ None => ElementKind::Quad4,
+ Some(k) => k
+ .as_str()
+ .and_then(ElementKind::from_name)
+ .filter(|k| k.dimension() == 2)
+ .ok_or_else(|| EngineError::Params("`kind` must be quad4, quad8, tri3 or tri6".into()))?,
+ };
+ let opts = modelcraft_mesh::quad::QuadOptions {
+ max_elements: optional_count(p, "max_elements", 1_000_000, 1..=PLANAR_MAX_ELEMENTS)?,
+ ..Default::default()
+ };
+ // The region's loops of curves, the body name, the face its elements fill, and each curve
+ // id's edge name.
+ let (loops, body, face, names): (Vec>, String, String, BTreeMap) = match (p.get("polygon"), p.get("body")) {
+ (Some(outer), None) => {
+ let mut rings = vec![polygon(outer, "polygon")?];
+ if let Some(h) = p.get("holes") {
+ let list = h
+ .as_array()
+ .filter(|a| a.len() <= PLANAR_MAX_HOLES)
+ .ok_or_else(|| EngineError::Params(format!("`holes` must be a list of at most {PLANAR_MAX_HOLES} polygons")))?;
+ for hole in list {
+ rings.push(polygon(hole, "holes")?);
+ }
+ }
+ let name = match p.get("name") {
+ None => "Plane".to_string(),
+ Some(n) => n.as_str().map(str::to_string).ok_or_else(|| EngineError::Params("`name` must be a body name".into()))?,
+ };
+ GeometryRef::body(name.clone())?;
+ let mut id = 0u32;
+ let mut names = BTreeMap::new();
+ let loops = rings
+ .into_iter()
+ .map(|ring| {
+ ring.into_iter()
+ .map(|pt| {
+ id += 1;
+ names.insert(id, format!("e{id}"));
+ Curve { id: Some(id), points: vec![pt] }
+ })
+ .collect()
+ })
+ .collect();
+ (loops, name, "area".to_string(), names)
+ }
+ (None, Some(_)) => {
+ let face_name = p.get("face").and_then(Value::as_str).ok_or_else(|| EngineError::Params("`face` must be a face name".into()))?;
+ let (shape, mesh) = super::geometry::body_mesh(s, p)?;
+ let face = shape.face_index(face_name).ok_or_else(|| EngineError::Params(format!("body `{}` has no face `{face_name}`", shape.name)))?;
+ let b = face_boundary(&mesh.positions, &mesh.triangles, &mesh.tri_face, face, mesh.weld_distance)?;
+ let names = shape.faces().iter().enumerate().map(|(i, f)| (i as u32, f.name.clone())).collect();
+ (b.loops, shape.name.clone(), face_name.to_string(), names)
+ }
+ _ => return Err(EngineError::Params("give either `polygon` (and `holes`) or `body` with `face` and `tolerance`".into())),
+ };
+ let m = mesh_region(&loops, kind, size, &opts)?;
+
+ let (first_node, first_element) = next_labels(s, m.nodes.len(), m.elements.len())?;
+ for (i, xy) in m.nodes.iter().enumerate() {
+ s.model.set_node(first_node + i as u32, [xy[0], xy[1], 0.0])?;
+ }
+ let region = s.model.geometry_tag(&GeometryRef::face(body.clone(), face.clone())?)?;
+ for (i, nodes) in m.elements.iter().enumerate() {
+ let label = first_element + i as u32;
+ s.model.set_element(label, kind, nodes.iter().map(|v| first_node + v).collect())?;
+ s.model.tag_element(label, region)?;
+ }
+ let mut tags: BTreeMap = BTreeMap::new();
+ for (e, edge, id) in &m.edges {
+ let Some(name) = names.get(id) else { continue };
+ let tag = match tags.get(id) {
+ Some(t) => *t,
+ None => {
+ let t = s.model.geometry_tag(&GeometryRef::new(body.clone(), Entity::Edge(name.clone()))?)?;
+ tags.insert(*id, t);
+ t
+ }
+ };
+ s.model.tag_element_edge(first_element + *e as u32, *edge, tag)?;
+ }
+ let polygon_input = p.get("polygon").is_some();
+ if polygon_input {
+ for (node, id) in &m.corners {
+ let t = s.model.geometry_tag(&GeometryRef::new(body.clone(), Entity::Vertex(format!("v{id}")))?)?;
+ s.model.tag_node(first_node + node, t)?;
+ }
+ }
+ let edges: Vec<&String> = tags.keys().filter_map(|id| names.get(id)).collect();
+ Ok(json!({
+ "body": body,
+ "face": face,
+ "kind": kind.name(),
+ "nodes": m.nodes.len(),
+ "elements": m.elements.len(),
+ "first_node": first_node,
+ "first_element": first_element,
+ "edges": edges,
+ }))
+}
+
fn node_add(s: &mut Session, p: &Value) -> Result {
let l = label(p, "label")?;
s.model.set_node(l, xyz(p, "xyz")?)?;
@@ -527,4 +670,80 @@ mod tests {
assert!(s.execute("mesh.tet.generate", &json!({"surface": 1})).is_err(), "an open surface bounds no volume");
assert_eq!(s.model.elements().len(), 1);
}
+
+ #[test]
+ fn planar_meshes_a_polygon_that_solves_as_a_2d_model() {
+ let mut s = Session::default();
+ let hole: Vec<[f64; 2]> = (0..12).map(|k| std::f64::consts::TAU * k as f64 / 12.0).map(|t| [5.0 + t.cos(), 2.0 + t.sin()]).collect();
+ let r = s
+ .execute(
+ "mesh.planar",
+ &json!({"polygon": [[0, 0], [10, 0], [10, 4], [0, 4]], "holes": [hole], "size": 0.5, "kind": "quad8", "name": "Strip"}),
+ )
+ .unwrap();
+ assert_eq!(r["kind"], "quad8");
+ let elements = r["elements"].as_u64().unwrap() as usize;
+ assert!(elements > 20, "{r}");
+ assert!(s.model.elements().values().all(|e| e.kind == modelcraft_model::ElementKind::Quad8));
+ assert!(s.model.nodes().values().all(|n| n.xyz[2] == 0.0));
+ // Sets on the tagged edges carry the boundary condition and the load.
+ for (cmd, params) in [
+ ("model.set.create", json!({"name": "all", "geometry": [{"body": "Strip", "face": "area"}]})),
+ ("model.set.create", json!({"name": "left", "geometry": [{"body": "Strip", "edge": "e4"}]})),
+ ("model.set.create", json!({"name": "right", "geometry": [{"body": "Strip", "edge": "e2"}]})),
+ ("model.material.create", json!({"name": "Steel", "elastic": {"type": "isotropic", "young": 210000.0, "poisson": 0.3}})),
+ ("model.section.create", json!({"name": "plate", "set": "all", "material": "Steel", "kind": {"type": "plane_stress", "thickness": 2.0}})),
+ ("model.boundary_condition.create", json!({"name": "hold", "type": "fixed", "target": {"set": "left"}})),
+ ("model.load.create", json!({"name": "pull", "type": "traction", "set": "right", "vector": [10.0, 0.0, 0.0]})),
+ ("model.step.create", json!({"name": "Pull", "procedure": {"type": "static"}, "boundary_conditions": ["hold"], "loads": ["pull"]})),
+ ] {
+ s.execute(cmd, ¶ms).unwrap_or_else(|e| panic!("{cmd}: {e}"));
+ }
+ let r = s.execute("solve.static", &json!({})).unwrap();
+ let total: f64 = r["reaction_total"][0].as_f64().unwrap();
+ // The traction over the right edge: 10 × 4 × 2.
+ assert!((total + 80.0).abs() < 1e-9, "{r}");
+ // A second mesh appends; failures change nothing.
+ let journal = s.journal.len();
+ let model = s.model.clone();
+ for bad in [
+ json!({"polygon": [[0, 0], [1, 0], [0, 1]], "size": 0}),
+ json!({"polygon": [[0, 0], [1, 0], [0, 1]], "size": 1e-9}),
+ json!({"polygon": [[0, 0], [1, 0], [0, 1]], "size": 0.1, "kind": "hex8"}),
+ json!({"polygon": [[0, 0], [1, 0]], "size": 0.1}),
+ json!({"polygon": [[0, 0], [1, 0], [1, 1], [0, 1]], "holes": [[[2, 2], [3, 2], [3, 3]]], "size": 0.1}),
+ json!({"polygon": [[0, 0], [1, 1], [1, 0], [0, 1]], "size": 0.1}),
+ json!({"polygon": [[0, 0], [1, 0], [0, 1]], "size": 0.1, "max_elements": 0}),
+ json!({"polygon": [[0, 0], [1, 0], [0, 1]], "size": 0.1, "name": ""}),
+ json!({"polygon": [[0, 0], [1, 0], [0, 1]], "size": 0.1, "body": "Body1", "face": "x"}),
+ json!({"body": "Body1", "face": "x", "tolerance": 0.1, "size": 1}),
+ json!({"size": 1}),
+ ] {
+ assert!(s.execute("mesh.planar", &bad).is_err(), "{bad}");
+ }
+ assert_eq!((s.journal.len(), &s.model), (journal, &model), "failures change nothing");
+ }
+
+ #[test]
+ fn planar_meshes_a_flat_face_of_a_design_with_edges_named_by_their_neighbours() {
+ use std::sync::Arc;
+
+ use modelcraft_geometry::open_design;
+ use modelcraft_geometry::sample::{Plate, plate_design};
+
+ let design = open_design(&plate_design(&Plate::default()).unwrap()).unwrap();
+ let names: Vec = design.bodies[0].faces().iter().map(|f| f.name.clone()).collect();
+ let mut s = Session { design: Some(Arc::new(design)), ..Session::default() };
+ let mut meshed = 0;
+ for face in &names {
+ let Ok(r) = s.execute("mesh.planar", &json!({"body": "Body1", "face": face, "tolerance": 0.2, "size": 4.0, "kind": "tri6"})) else {
+ continue;
+ };
+ meshed += 1;
+ let edges: Vec<&str> = r["edges"].as_array().unwrap().iter().map(|e| e.as_str().unwrap()).collect();
+ // The plate's top and bottom faces border its four sides and the hole.
+ assert!(edges.len() >= 5 && edges.iter().all(|e| names.iter().any(|n| n == e)), "{face}: {edges:?}");
+ }
+ assert_eq!(meshed, 2, "the top and bottom faces are flat and parallel to XY; the sides and hole are not");
+ }
}
diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs
index 8ed73e4..7cf562d 100644
--- a/crates/engine/src/cmd/mod.rs
+++ b/crates/engine/src/cmd/mod.rs
@@ -30,6 +30,7 @@ pub static ALL: &[CommandSpec] = &[
mesh::QUALITY,
mesh::DELAUNAY_POINTS,
mesh::SURFACE_CHECK,
+ mesh::PLANAR,
mesh::TET_GENERATE,
tags::TAG,
tags::TAG_CLEAR,
diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs
index f7a829e..d47392d 100644
--- a/crates/engine/src/lib.rs
+++ b/crates/engine/src/lib.rs
@@ -330,6 +330,13 @@ mod tests {
json!({"surface": -1}),
json!({"surface": 1e300}),
json!({"surface": "0", "weld": -1, "feature_angle": 1e300}),
+ json!({"polygon": vec![[0.0, 0.0]; 200_000], "size": 1}),
+ json!({"polygon": [[0, 0], [1e308, 0], [0, 1e308]], "size": 1e-300}),
+ json!({"polygon": [[0, 0], [1, 0], [0, 1]], "holes": vec![[[0.1, 0.1], [0.2, 0.1], [0.1, 0.2]]; 20_000], "size": 0.5}),
+ json!({"polygon": [[0, 0], [1, 0], [0, 1]], "holes": [[[0.1, 0.1], [0.2, 0.1], [0.1, 0.2]], [[0.1, 0.1], [0.2, 0.1], [0.1, 0.2]]], "size": 0.2}),
+ json!({"polygon": [["a", 0], [1, null], {}], "size": "big", "kind": 7}),
+ json!({"polygon": [[0, 0], [1, 0], [0, 1]], "size": 1e300, "kind": "tri6", "max_elements": 1e300}),
+ json!({"body": "Body1", "face": "", "tolerance": 0.05, "size": 1}),
// Model definitions.
json!({"name": "x".repeat(100_000)}),
json!({"name": "A\u{0}"}),
@@ -345,6 +352,12 @@ mod tests {
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", "set": "all", "material": "Steel", "kind": {"type": "plane_stress", "thickness": -1}}),
+ json!({"name": "A", "set": "all", "material": "Steel", "kind": {"type": "plane_strain", "thickness": 1e308}}),
+ json!({"name": "A", "set": "all", "material": "Steel", "kind": {"type": "axisymmetric", "axis": "x"}}),
+ json!({"name": "A", "type": "centrifugal", "angular_velocity": 1e300, "point": [0, 0, 0], "axis": [0, 1, 0]}),
+ json!({"name": "A", "type": "centrifugal", "angular_velocity": 1, "point": [0, 0, 0], "axis": [0, 0, 0]}),
+ json!({"name": "A", "edges": [{"element": u32::MAX, "edge": 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]}),
diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs
index d61dc9b..51c3f70 100644
--- a/crates/explicit/src/system.rs
+++ b/crates/explicit/src/system.rs
@@ -368,6 +368,17 @@ impl System {
}
});
}
+ Load::Centrifugal { angular_velocity, point, axis, set } => {
+ let chosen: Option<&[u32]> = set.as_ref().map(|s| resolved.sets.get(s).map(|r| r.elements.as_slice()).unwrap_or_default());
+ let xyz = &self.xyz;
+ self.each_element(chosen, |nodes, weights, rho| {
+ for (node, w) in nodes.iter().zip(weights) {
+ let x = xyz.get(*node).copied().unwrap_or_default();
+ let a = modelcraft_model::centrifugal_acceleration(*angular_velocity, *point, *axis, x);
+ add_at(&mut f, *node, &a.map(|c| rho * w * c));
+ }
+ });
+ }
}
}
self.external = f;
diff --git a/crates/fem/src/element.rs b/crates/fem/src/element.rs
index 2ef905a..9b3bf4b 100644
--- a/crates/fem/src/element.rs
+++ b/crates/fem/src/element.rs
@@ -34,6 +34,7 @@
use modelcraft_model::ElementKind;
+use crate::basis::Basis;
use crate::dense::{Cholesky, Matrix, inverse3, positive_inverse};
use crate::dof::{DofSet, ElementDofs, Internal};
use crate::material::{Elasticity, stress};
@@ -47,6 +48,9 @@ pub enum Family {
Tetrahedron,
Wedge,
Hexahedron,
+ /// The 2D elements of plane and axisymmetric models ([`crate::plane`]).
+ Triangle,
+ Quadrilateral,
}
/// One row of the element table: everything the kernels need to know about an element kind.
@@ -386,15 +390,7 @@ impl<'a> Solid<'a> {
/// The lumped (diagonal) mass per node: the consistent diagonal scaled to the element's mass.
pub fn lumped_mass(&self, density: f64) -> Result> {
- let m = self.mass(density)?;
- // The consistent entries sum to the element's mass (the shape functions sum to one).
- let total: f64 = m.data.iter().sum();
- let diag: Vec = (0..m.n).map(|i| m.get(i, i)).collect();
- let sum: f64 = diag.iter().sum();
- if !(sum > 0.0) {
- return Err(self.degenerate());
- }
- Ok(diag.iter().map(|d| d * total / sum).collect())
+ lump(&self.mass(density)?).ok_or_else(|| self.degenerate())
}
/// Consistent nodal forces of a uniform force per volume: fᵢ = ∫ Nᵢ b dV.
@@ -411,6 +407,22 @@ impl<'a> Solid<'a> {
Ok(f)
}
+ /// Consistent nodal forces of a force per volume that varies with position (such as a
+ /// centrifugal load): fᵢ = ∫ Nᵢ b(x) dV.
+ pub fn body_force_field(&self, b: impl Fn([f64; 3]) -> [f64; 3]) -> Result> {
+ let mut f = vec![[0.0; 3]; self.ty.nodes];
+ for p in self.ty.mass.points() {
+ let g = self.gradients(p.xi)?;
+ let b = b(self.position(p.xi));
+ for (fi, n) in f.iter_mut().zip(g.shape.values()) {
+ for k in 0..3 {
+ fi[k] += n * b[k] * g.det * p.weight;
+ }
+ }
+ }
+ Ok(f)
+ }
+
/// Strain and stress at every integration point for the nodal displacements `u`.
pub fn stresses(&self, d: &Elasticity, u: &[[f64; 3]]) -> Result> {
if u.len() != self.ty.nodes {
@@ -475,7 +487,16 @@ impl<'a> Solid<'a> {
}
/// Make a matrix exactly symmetric (the mean of it and its transpose), so assembly can rely on it.
-fn symmetrize(k: &mut Matrix) {
+/// HRZ lumping of a consistent mass matrix: its diagonal scaled to the element's mass (the
+/// entries sum to it, the shape functions summing to one); `None` when the diagonal has no mass.
+pub(crate) fn lump(m: &Matrix) -> Option> {
+ let total: f64 = m.data.iter().sum();
+ let diag: Vec = (0..m.n).map(|i| m.get(i, i)).collect();
+ let sum: f64 = diag.iter().sum();
+ (sum > 0.0).then(|| diag.iter().map(|d| d * total / sum).collect())
+}
+
+pub(crate) fn symmetrize(k: &mut Matrix) {
let n = k.n;
for i in 0..n {
for j in i + 1..n {
@@ -589,6 +610,8 @@ pub fn extrapolation(ty: &ElementType, rule: Rule) -> Result {
Family::Tetrahedron => shape::tet4(xi),
Family::Wedge => shape::wedge6(xi),
Family::Hexahedron => shape::hex8(xi),
+ Family::Triangle => Basis::Tri3.eval(xi),
+ Family::Quadrilateral => Basis::Quad4.eval(xi),
};
let mut corner_rows = vec![0.0; c * g];
if g < c {
diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs
index 53e0ec8..b22806d 100644
--- a/crates/fem/src/lib.rs
+++ b/crates/fem/src/lib.rs
@@ -15,6 +15,8 @@
//! - [`material`]: isotropic and orthotropic elasticity matrices, von Mises stress;
//! - [`element`]: the element table, stiffness (full, reduced, incompatible modes), consistent
//! and lumped mass, body forces, integration-point stresses and their extrapolation to nodes;
+//! - [`plane`]: 2D elements (plane stress, plane strain, axisymmetric) on tri3, tri6, quad4 and
+//! quad8, with two displacement components per node;
//! - [`face`]: consistent nodal forces of pressure on element faces;
//! - [`uniform`]: uniform-strain (one-point) hexahedron and tetrahedron kinematics, hourglass
//! shape vectors and element frequency bounds for explicit dynamics;
@@ -35,6 +37,7 @@ pub mod element;
pub mod face;
pub mod jet;
pub mod material;
+pub mod plane;
pub mod quadrature;
pub mod shape;
pub mod uniform;
diff --git a/crates/fem/src/plane.rs b/crates/fem/src/plane.rs
new file mode 100644
index 0000000..e925b2b
--- /dev/null
+++ b/crates/fem/src/plane.rs
@@ -0,0 +1,701 @@
+//! 2D continuum elements: plane stress, plane strain and axisymmetric.
+//!
+//! The elements lie in the XY plane and have two displacement components per node, u (along X)
+//! and v (along Y). The strain is the 4-vector (εxx, εyy, εzz, γxy), where the third component is
+//! the out-of-plane one:
+//! - **plane strain** (a slice of a long body): εzz = 0, and σzz = D₃₁εxx + D₃₂εyy + D₃₄γxy is
+//! reported;
+//! - **plane stress** (a thin plate loaded in its plane): σzz = 0, the elasticity is condensed to
+//! D̃ᵢⱼ = Dᵢⱼ − Dᵢ₃D₃ⱼ/D₃₃ and εzz = −(D₃₁εxx + D₃₂εyy + D₃₄γxy)/D₃₃ is reported;
+//! - **axisymmetric** (a body of revolution about the Y axis): x is the radius r, y the axial
+//! coordinate, and εzz is the hoop strain εθθ = u/r; volumes are 2πr dA, so forces, reactions
+//! and masses are totals over the full circumference.
+//!
+//! D is the leading 4 × 4 block (xx, yy, zz, xy) of the 3D elasticity, so orthotropic materials
+//! with axes along the global axes work as they do in solids; stresses and strains are reported
+//! as the solid elements' 6-vectors (xx, yy, zz, xy, xz, yz) with xz = yz = 0 (axisymmetric: rr,
+//! axial, hoop, r-axial). Plane elements integrate over their thickness t, so stiffness, mass and
+//! loads are those of the whole thickness. (O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, *The
+//! Finite Element Method: Its Basis and Fundamentals*, 7th ed., Butterworth-Heinemann 2013,
+//! ch. 6 (plane stress and strain) and §6.6 (axisymmetric); T. J. R. Hughes, *The Finite Element
+//! Method*, Dover 2000, §2.7 and §3.5.)
+//!
+//! Elements and rules: tri3 (one point), tri6 (three points, degree 2), quad4 (2 × 2) and quad8
+//! (3 × 3, or 2 × 2 reduced). Reduced quad4 needs hourglass control and is not offered. Quad4 has
+//! an **incompatible modes** formulation (QM6): the internal modes (1 − ξ²) and (1 − η²) in u
+//! and v (E. L. Wilson, R. L. Taylor, W. P. Doherty, J. Ghaboussi, *Incompatible displacement
+//! models*, 1973), their natural derivatives mapped with the centre Jacobian and scaled by
+//! det J₀/det J (R. L. Taylor, P. J. Beresford, E. L. Wilson, *A non-conforming element for
+//! stress analysis*, IJNME 10(6):1211–1219, 1976), condensed out of the element. In an
+//! axisymmetric element the modes also strain the hoop direction (P/r), and Taylor's scaling no
+//! longer makes their strains integrate to zero over the element's volume; their strain
+//! operator is therefore shifted by its volume mean, the patch-test correction of E. L. Wilson
+//! and A. Ibrahimbegovic (*Use of incompatible displacement modes for the calculation of element
+//! stiffnesses or stresses*, Finite Elements in Analysis and Design 7(3):229–241, 1990), which
+//! leaves the plane case unchanged (its mean is zero) and makes the axisymmetric element pass the
+//! patch test.
+//!
+//! On the axis (r = 0) the hoop strain u/r is evaluated only at integration points, which lie
+//! inside the element (r > 0) for every element with some node off the axis; the solver holds
+//! the radial displacement of nodes on the axis at zero, as symmetry requires.
+
+use modelcraft_model::{ElementKind, Extent, SectionKind};
+
+use crate::basis::Basis;
+use crate::dense::{Cholesky, Matrix};
+use crate::dof::{DofSet, ElementDofs, Internal};
+use crate::element::{ElementType, Family, Formulation, PointStress, extrapolation, lump, symmetrize};
+use crate::material::Elasticity;
+use crate::quadrature::{Point, Rule};
+use crate::shape::{MAX_NODES, Shape};
+use crate::{Error, Result};
+
+/// Which 2D idealisation an element follows.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub enum Theory {
+ PlaneStress {
+ thickness: f64,
+ },
+ PlaneStrain {
+ thickness: f64,
+ },
+ /// About the Y axis; x is the radius.
+ Axisymmetric,
+}
+
+/// The 4 × 4 elasticity of a 2D element, in the order (xx, yy, zz, xy).
+pub type Elasticity2 = [[f64; 4]; 4];
+
+impl Theory {
+ /// The theory of a 2D section (`None` for solid, shell, beam and rod sections).
+ pub fn from_section(kind: &SectionKind) -> Option {
+ match *kind {
+ SectionKind::PlaneStress { thickness } => Some(Theory::PlaneStress { thickness }),
+ SectionKind::PlaneStrain { thickness } => Some(Theory::PlaneStrain { thickness }),
+ SectionKind::Axisymmetric {} => Some(Theory::Axisymmetric),
+ _ => None,
+ }
+ }
+
+ /// How the element extends out of its plane (thickness, or revolved).
+ pub fn extent(&self) -> Extent {
+ match *self {
+ Theory::PlaneStress { thickness } | Theory::PlaneStrain { thickness } => Extent::Thickness(thickness),
+ Theory::Axisymmetric => Extent::Revolved,
+ }
+ }
+
+ /// The 4 × 4 elasticity of this theory from the 3D one.
+ pub fn elasticity(&self, d: &Elasticity) -> Result {
+ let mut out = [[0.0; 4]; 4];
+ for i in 0..4 {
+ for j in 0..4 {
+ out[i][j] = d[i][j];
+ }
+ }
+ if let Theory::PlaneStress { .. } = self {
+ let d33 = d[2][2];
+ if !(d33 > 0.0) {
+ return Err(Error::Invalid("the material has no stiffness through the thickness".into()));
+ }
+ for i in 0..4 {
+ for j in 0..4 {
+ out[i][j] = if i == 2 || j == 2 { 0.0 } else { d[i][j] - d[i][2] * d[2][j] / d33 };
+ }
+ }
+ }
+ Ok(out)
+ }
+}
+
+/// The 2D element table: chosen by the section kind (a 2D section on tri3, tri6, quad4 or quad8),
+/// separately from shells on the same element kinds.
+pub static PLANE_TYPES: [ElementType; 4] = [
+ ElementType {
+ kind: ElementKind::Tri3,
+ family: Family::Triangle,
+ nodes: 3,
+ corners: 3,
+ shape: |xi| Basis::Tri3.eval(xi),
+ full: Rule::TriCentroid,
+ reduced: None,
+ mass: Rule::TriConical(3),
+ incompatible_modes: false,
+ },
+ ElementType {
+ kind: ElementKind::Tri6,
+ family: Family::Triangle,
+ nodes: 6,
+ corners: 3,
+ shape: |xi| Basis::Tri6.eval(xi),
+ full: Rule::Tri3,
+ reduced: None,
+ mass: Rule::TriConical(4),
+ incompatible_modes: false,
+ },
+ ElementType {
+ kind: ElementKind::Quad4,
+ family: Family::Quadrilateral,
+ nodes: 4,
+ corners: 4,
+ shape: |xi| Basis::Quad4.eval(xi),
+ full: Rule::Quad(2),
+ reduced: None,
+ mass: Rule::Quad(2),
+ incompatible_modes: true,
+ },
+ ElementType {
+ kind: ElementKind::Quad8,
+ family: Family::Quadrilateral,
+ nodes: 8,
+ corners: 4,
+ shape: |xi| Basis::Quad8.eval(xi),
+ full: Rule::Quad(3),
+ reduced: Some(Rule::Quad(2)),
+ mass: Rule::Quad(3),
+ incompatible_modes: false,
+ },
+];
+
+/// The 2D table row of an element kind.
+pub fn plane_type(kind: ElementKind) -> Result<&'static ElementType> {
+ PLANE_TYPES
+ .iter()
+ .find(|t| t.kind == kind)
+ .ok_or_else(|| Error::NotSupported(format!("{} elements in a 2D section (2D sections take tri3, tri6, quad4 and quad8)", kind.name())))
+}
+
+/// The formulation a 2D section asks for, if this element type has it.
+pub fn formulation(kind: ElementKind, f: &modelcraft_model::Formulation) -> Result {
+ let ty = plane_type(kind)?;
+ if f.hybrid {
+ return Err(Error::NotSupported("hybrid (mixed displacement–pressure) 2D elements".into()));
+ }
+ if f.incompatible_modes {
+ return if ty.incompatible_modes {
+ Ok(Formulation::IncompatibleModes)
+ } else {
+ Err(Error::Invalid(format!("incompatible modes apply to quad4 elements, not {}", kind.name())))
+ };
+ }
+ match f.integration {
+ modelcraft_model::Integration::Full => Ok(Formulation::Full),
+ modelcraft_model::Integration::Reduced if ty.reduced.is_some() => Ok(Formulation::Reduced),
+ modelcraft_model::Integration::Reduced => Err(Error::NotSupported(format!(
+ "reduced integration of 2D {} elements (quad4 needs hourglass control; use full integration or incompatible modes)",
+ kind.name()
+ ))),
+ }
+}
+
+/// One column of the strain operator: (∂N/∂x, ∂N/∂y, hoop N/r), for a node or an internal mode.
+type Column = [f64; 3];
+
+/// What the kernels need at one integration point.
+struct At {
+ shape: Shape,
+ /// Strain columns of the nodes.
+ g: [Column; MAX_NODES],
+ det: f64,
+ x: [f64; 3],
+ /// det J times the point's weight times the out-of-plane extent: the point's share of the
+ /// volume.
+ dv: f64,
+}
+
+/// The incompatible modes' count of internal degrees of freedom: two modes, two directions.
+const INCOMPATIBLE_DOFS: usize = 4;
+
+/// A 2D element: its type, node coordinates (z = 0), formulation and theory.
+#[derive(Debug, Clone, Copy)]
+pub struct Plane<'a> {
+ pub ty: &'static ElementType,
+ xyz: &'a [[f64; 3]],
+ formulation: Formulation,
+ theory: Theory,
+}
+
+/// out[2i + r][2j + c] += w (Bᵢᵀ D Bⱼ)[r][c] over every pair of columns (row-major, `cols` wide),
+/// computed for j ≥ i and mirrored, so the result is exactly symmetric.
+fn accumulate(out: &mut [f64], columns: &[Column], d: &Elasticity2, w: f64) {
+ let cols = 2 * columns.len();
+ let block = |g: &Column| -> [[f64; 2]; 4] { [[g[0], 0.0], [0.0, g[1]], [g[2], 0.0], [g[1], g[0]]] };
+ for (j, gj) in columns.iter().enumerate() {
+ let bj = block(gj);
+ let mut db = [[0.0; 2]; 4];
+ for s in 0..4 {
+ for c in 0..2 {
+ db[s][c] = (0..4).map(|t| d[s][t] * bj[t][c]).sum();
+ }
+ }
+ for (i, gi) in columns.iter().enumerate().take(j + 1) {
+ let bi = block(gi);
+ for r in 0..2 {
+ for c in if i == j { r } else { 0 }..2 {
+ let v = w * (0..4).map(|s| bi[s][r] * db[s][c]).sum::();
+ let (row, col) = (2 * i + r, 2 * j + c);
+ if let Some(x) = out.get_mut(row * cols + col) {
+ *x += v;
+ }
+ if row != col
+ && let Some(x) = out.get_mut(col * cols + row)
+ {
+ *x += v;
+ }
+ }
+ }
+ }
+ }
+}
+
+impl<'a> Plane<'a> {
+ pub fn new(kind: ElementKind, xyz: &'a [[f64; 3]], formulation: Formulation, theory: Theory) -> Result> {
+ let ty = plane_type(kind)?;
+ if xyz.len() != ty.nodes {
+ return Err(Error::Invalid(format!("a {} element has {} nodes, got {}", kind.name(), ty.nodes, xyz.len())));
+ }
+ modelcraft_model::check_coordinates(xyz.as_flattened()).map_err(|e| Error::Invalid(e.to_string()))?;
+ if xyz.iter().any(|p| p[2] != 0.0) {
+ return Err(Error::Invalid(format!("a 2D {} element must lie in the XY plane (z = 0)", kind.name())));
+ }
+ if theory == Theory::Axisymmetric && xyz.iter().any(|p| p[0] < 0.0) {
+ return Err(Error::Invalid("an axisymmetric element must lie in x ≥ 0 (x is the radius)".into()));
+ }
+ if let Theory::PlaneStress { thickness } | Theory::PlaneStrain { thickness } = theory
+ && !(thickness > 0.0 && thickness.is_finite())
+ {
+ return Err(Error::Invalid(format!("the thickness must be a positive number, got {thickness}")));
+ }
+ let ok = match formulation {
+ Formulation::Full => true,
+ Formulation::Reduced => ty.reduced.is_some(),
+ Formulation::IncompatibleModes => ty.incompatible_modes,
+ };
+ if !ok {
+ return Err(Error::NotSupported(format!("{formulation:?} formulation of 2D {} elements", kind.name())));
+ }
+ Ok(Plane { ty, xyz, formulation, theory })
+ }
+
+ /// The element's unknowns: u and v at every node ([`DofSet::PLANE`]), and the incompatible
+ /// modes it condenses out.
+ pub fn dofs(&self) -> ElementDofs {
+ let incompatible_modes = if self.formulation == Formulation::IncompatibleModes { INCOMPATIBLE_DOFS } else { 0 };
+ ElementDofs { node: DofSet::PLANE, internal: Internal { incompatible_modes, pressure: 0 } }
+ }
+
+ /// The stiffness integration rule.
+ pub fn rule(&self) -> Rule {
+ match (self.formulation, self.ty.reduced) {
+ (Formulation::Reduced, Some(r)) => r,
+ _ => self.ty.full,
+ }
+ }
+
+ fn degenerate(&self) -> Error {
+ Error::Degenerate(format!("a 2D {} element is inverted, clockwise or degenerate (its Jacobian is not positive)", self.ty.kind.name()))
+ }
+
+ /// The 2 × 2 Jacobian J[a][b] = ∂x_b/∂ξ_a and its determinant.
+ fn jacobian(&self, s: &Shape) -> ([[f64; 2]; 2], f64) {
+ let mut j = [[0.0; 2]; 2];
+ for (d, x) in s.derivatives().iter().zip(self.xyz) {
+ for a in 0..2 {
+ for b in 0..2 {
+ j[a][b] += d[a] * x[b];
+ }
+ }
+ }
+ (j, j[0][0] * j[1][1] - j[0][1] * j[1][0])
+ }
+
+ fn at(&self, p: &Point) -> Result {
+ let basis = Basis::of(self.ty.kind);
+ let shape = basis.eval(p.xi);
+ let grads = basis.gradients(&shape, self.xyz).map_err(|_| self.degenerate())?;
+ let det = grads.det;
+ let mut x = [0.0; 3];
+ for (n, xi) in shape.values().iter().zip(self.xyz) {
+ for k in 0..2 {
+ x[k] += n * xi[k];
+ }
+ }
+ let axisymmetric = self.theory == Theory::Axisymmetric;
+ if axisymmetric && !(x[0] > 0.0) {
+ return Err(Error::Degenerate(format!("an axisymmetric {} element has an integration point on the axis", self.ty.kind.name())));
+ }
+ let mut g = [[0.0; 3]; MAX_NODES];
+ for ((gi, d), n) in g.iter_mut().zip(grads.gradients()).zip(shape.values()) {
+ *gi = [d[0], d[1], if axisymmetric { n / x[0] } else { 0.0 }];
+ }
+ let dv = det * p.weight * self.theory.extent().at(x);
+ Ok(At { shape, g, det, x, dv })
+ }
+
+ /// The element's volume (area times thickness, or the volume of the ring).
+ pub fn volume(&self) -> Result {
+ let mut v = 0.0;
+ for p in self.ty.mass.points() {
+ v += self.at(&p)?.dv;
+ }
+ Ok(v)
+ }
+
+ /// The incompatible modes' strain columns at every stiffness point: Taylor's mapped
+ /// derivatives (and the hoop value P/r when axisymmetric), shifted by their volume mean.
+ fn incompatible_columns(&self, points: &[(Point, At)]) -> Result> {
+ let s0 = (self.ty.shape)([0.0; 3]);
+ let (j0, det0) = self.jacobian(&s0);
+ if !(det0 > 0.0 && det0.is_finite()) {
+ return Err(self.degenerate());
+ }
+ let axisymmetric = self.theory == Theory::Axisymmetric;
+ let mut out = Vec::with_capacity(points.len());
+ let mut mean = [[0.0; 3]; 2];
+ let mut volume = 0.0;
+ for (p, a) in points {
+ let mut cols = [[0.0; 3]; 2];
+ for (k, col) in cols.iter_mut().enumerate() {
+ // P_k = 1 − ξ_k²: ∂P/∂ξ_k = −2ξ_k, the other derivative zero.
+ let mut dp = [0.0; 2];
+ dp[k] = -2.0 * p.xi[k];
+ // (det J₀ / det J) J₀⁻¹ ∂P/∂ξ, with J₀⁻¹ = adj(J₀) / det J₀.
+ col[0] = (j0[1][1] * dp[0] - j0[0][1] * dp[1]) / a.det;
+ col[1] = (-j0[1][0] * dp[0] + j0[0][0] * dp[1]) / a.det;
+ col[2] = if axisymmetric { (1.0 - p.xi[k] * p.xi[k]) / a.x[0] } else { 0.0 };
+ for c in 0..3 {
+ mean[k][c] += col[c] * a.dv;
+ }
+ }
+ volume += a.dv;
+ out.push(cols);
+ }
+ if !(volume > 0.0) {
+ return Err(self.degenerate());
+ }
+ for cols in &mut out {
+ for k in 0..2 {
+ for c in 0..3 {
+ cols[k][c] -= mean[k][c] / volume;
+ }
+ }
+ }
+ Ok(out)
+ }
+
+ /// The integration points with their data.
+ fn points(&self) -> Result> {
+ self.rule().points().into_iter().map(|p| self.at(&p).map(|a| (p, a))).collect()
+ }
+
+ /// The stiffness over nodes and, with incompatible modes, the internal modes after them:
+ /// (2n + 4) × (2n + 4).
+ fn full_stiffness(&self, d: &Elasticity2) -> Result<(Matrix, Vec<(Point, At)>, Option>)> {
+ let n = self.ty.nodes;
+ let points = self.points()?;
+ let modes = if self.formulation == Formulation::IncompatibleModes { Some(self.incompatible_columns(&points)?) } else { None };
+ let width = 2 * n + modes.as_ref().map_or(0, |_| INCOMPATIBLE_DOFS);
+ let mut k = Matrix::zeros(width);
+ let mut columns: Vec = Vec::with_capacity(n + 2);
+ for (i, (_, a)) in points.iter().enumerate() {
+ columns.clear();
+ columns.extend_from_slice(a.g.get(..n).unwrap_or_default());
+ if let Some(m) = modes.as_ref().and_then(|m| m.get(i)) {
+ columns.extend_from_slice(m);
+ }
+ accumulate(&mut k.data, &columns, d, a.dv);
+ }
+ Ok((k, points, modes))
+ }
+
+ /// The stiffness matrix, 2n × 2n, degrees of freedom ordered node by node (u, v).
+ pub fn stiffness(&self, d: &Elasticity) -> Result {
+ let d2 = self.theory.elasticity(d)?;
+ let (full, _, modes) = self.full_stiffness(&d2)?;
+ if modes.is_none() {
+ return Ok(full);
+ }
+ // K = K_uu − K_ua K_aa⁻¹ K_au.
+ let n = 2 * self.ty.nodes;
+ let kaa = sub_block(&full, n, n, INCOMPATIBLE_DOFS);
+ let chol = Cholesky::new(&kaa)?;
+ let mut k = sub_block(&full, 0, 0, n);
+ for i in 0..n {
+ let mut col: Vec = (0..INCOMPATIBLE_DOFS).map(|a| full.get(n + a, i)).collect();
+ chol.solve(&mut col);
+ for j in 0..n {
+ let v: f64 = (0..INCOMPATIBLE_DOFS).map(|a| full.get(j, n + a) * col.get(a).copied().unwrap_or(0.0)).sum();
+ k.add(j, i, -v);
+ }
+ }
+ symmetrize(&mut k);
+ Ok(k)
+ }
+
+ /// The consistent mass matrix per displacement component, n × n: m_ij = ∫ ρ Nᵢ Nⱼ dV.
+ pub fn mass(&self, density: f64) -> Result {
+ let n = self.ty.nodes;
+ let mut m = Matrix::zeros(n);
+ for p in self.ty.mass.points() {
+ let a = self.at(&p)?;
+ let w = density * a.dv;
+ let v = a.shape.values();
+ for (i, x) in v.iter().enumerate() {
+ for (j, y) in v.iter().enumerate() {
+ m.add(i, j, w * x * y);
+ }
+ }
+ }
+ Ok(m)
+ }
+
+ /// The lumped (diagonal) mass per node: the consistent diagonal scaled to the element's mass
+ /// (Hinton, Rock and Zienkiewicz 1976).
+ pub fn lumped_mass(&self, density: f64) -> Result> {
+ lump(&self.mass(density)?).ok_or_else(|| self.degenerate())
+ }
+
+ /// Consistent nodal forces of a force per volume b(x) in the plane: fᵢ = ∫ Nᵢ b(x) dV (the Z
+ /// component of `b` is ignored).
+ pub fn body_force_field(&self, b: impl Fn([f64; 3]) -> [f64; 3]) -> Result> {
+ let mut f = vec![[0.0; 2]; self.ty.nodes];
+ for p in self.ty.mass.points() {
+ let a = self.at(&p)?;
+ let b = b(a.x);
+ for (fi, n) in f.iter_mut().zip(a.shape.values()) {
+ for k in 0..2 {
+ fi[k] += n * b[k] * a.dv;
+ }
+ }
+ }
+ Ok(f)
+ }
+
+ /// Strain and stress at every integration point for the nodal displacements `u` (u, v).
+ pub fn stresses(&self, d: &Elasticity, u: &[[f64; 2]]) -> Result> {
+ if u.len() != self.ty.nodes {
+ return Err(Error::Invalid(format!("{} displacements for a {}-node element", u.len(), self.ty.nodes)));
+ }
+ let d2 = self.theory.elasticity(d)?;
+ let (full, points, modes) =
+ if self.formulation == Formulation::IncompatibleModes { self.full_stiffness(&d2)? } else { (Matrix::zeros(0), self.points()?, None) };
+ // Incompatible modes: α = −K_aa⁻¹ K_au u.
+ let alpha = if modes.is_some() {
+ let n = 2 * self.ty.nodes;
+ let flat: Vec = u.iter().flatten().copied().collect();
+ let mut rhs: Vec =
+ (0..INCOMPATIBLE_DOFS).map(|a| -flat.iter().enumerate().map(|(i, x)| full.get(n + a, i) * x).sum::()).collect();
+ Cholesky::new(&sub_block(&full, n, n, INCOMPATIBLE_DOFS))?.solve(&mut rhs);
+ rhs
+ } else {
+ Vec::new()
+ };
+ let mut out = Vec::with_capacity(points.len());
+ for (i, (_, a)) in points.iter().enumerate() {
+ let mut e = [0.0; 4];
+ let mut add = |g: &Column, uv: [f64; 2]| {
+ e[0] += g[0] * uv[0];
+ e[1] += g[1] * uv[1];
+ e[2] += g[2] * uv[0];
+ e[3] += g[1] * uv[0] + g[0] * uv[1];
+ };
+ for (g, uv) in a.g.iter().zip(u) {
+ add(g, *uv);
+ }
+ if let Some(cols) = modes.as_ref().and_then(|m| m.get(i)) {
+ for (k, g) in cols.iter().enumerate() {
+ add(g, [alpha.get(2 * k).copied().unwrap_or(0.0), alpha.get(2 * k + 1).copied().unwrap_or(0.0)]);
+ }
+ }
+ if let Theory::PlaneStress { .. } = self.theory {
+ e[2] = -(d[2][0] * e[0] + d[2][1] * e[1] + d[2][3] * e[3]) / d[2][2];
+ }
+ // Plane stress: D̃ has no zz row or column, so the εzz just recovered adds nothing.
+ let s: [f64; 4] = std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * e[c]).sum());
+ out.push(PointStress { xyz: a.x, strain: [e[0], e[1], e[2], e[3], 0.0, 0.0], stress: [s[0], s[1], s[2], s[3], 0.0, 0.0] });
+ }
+ Ok(out)
+ }
+
+ /// Extrapolate values at the integration points to the element's nodes.
+ pub fn extrapolate(&self, values: &[[f64; 6]]) -> Result> {
+ let e = extrapolation(self.ty, self.rule())?;
+ if values.len() != e.points {
+ return Err(Error::Invalid(format!("{} values for {} integration points", values.len(), e.points)));
+ }
+ Ok(e.apply(values))
+ }
+}
+
+/// The `size` × `size` block of `m` starting at (`row`, `col`).
+fn sub_block(m: &Matrix, row: usize, col: usize, size: usize) -> Matrix {
+ let mut out = Matrix::zeros(size);
+ for i in 0..size {
+ for j in 0..size {
+ out.add(i, j, m.get(row + i, col + j));
+ }
+ }
+ out
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+ use crate::material::isotropic;
+
+ /// A distorted element of `kind` with corners in x ∈ [1, 3] (so axisymmetric elements stay off
+ /// the axis); mid-side nodes off the chord midpoints by a little, so quadratic elements are
+ /// curved.
+ fn element(kind: ElementKind) -> Vec<[f64; 3]> {
+ let corners: Vec<[f64; 3]> = match kind {
+ ElementKind::Tri3 | ElementKind::Tri6 => vec![[1.0, 0.2, 0.0], [3.0, 0.0, 0.0], [1.6, 1.9, 0.0]],
+ _ => vec![[1.0, 0.1, 0.0], [3.0, 0.0, 0.0], [2.7, 2.2, 0.0], [1.2, 1.6, 0.0]],
+ };
+ let mut xyz = corners.clone();
+ for (k, e) in kind.edges().iter().enumerate() {
+ if let [a, b, _] = **e {
+ let bump = 0.05 * (k as f64 + 1.0);
+ xyz.push([(corners[a][0] + corners[b][0]) / 2.0 + bump, (corners[a][1] + corners[b][1]) / 2.0 - bump, 0.0]);
+ }
+ }
+ xyz
+ }
+
+ fn variants() -> Vec<(ElementKind, Formulation)> {
+ vec![
+ (ElementKind::Tri3, Formulation::Full),
+ (ElementKind::Tri6, Formulation::Full),
+ (ElementKind::Quad4, Formulation::Full),
+ (ElementKind::Quad4, Formulation::IncompatibleModes),
+ (ElementKind::Quad8, Formulation::Full),
+ (ElementKind::Quad8, Formulation::Reduced),
+ ]
+ }
+
+ const THEORIES: [Theory; 3] = [Theory::PlaneStress { thickness: 0.7 }, Theory::PlaneStrain { thickness: 1.0 }, Theory::Axisymmetric];
+
+ #[test]
+ fn rigid_motion_strains_nothing_and_matrices_are_symmetric() {
+ let d = isotropic(1000.0, 0.3).unwrap();
+ for (kind, f) in variants() {
+ let xyz = element(kind);
+ for theory in THEORIES {
+ let e = Plane::new(kind, &xyz, f, theory).unwrap();
+ let k = e.stiffness(&d).unwrap();
+ let n = k.n;
+ for i in 0..n {
+ assert!(k.get(i, i) > 0.0);
+ for j in 0..n {
+ assert_eq!(k.get(i, j), k.get(j, i), "{kind:?} {f:?} {theory:?}");
+ }
+ }
+ // Translation along Y always; along X and the rotation only in plane models.
+ let mut motions: Vec> = vec![xyz.iter().flat_map(|_| [0.0, 1.0]).collect()];
+ if theory != Theory::Axisymmetric {
+ motions.push(xyz.iter().flat_map(|_| [1.0, 0.0]).collect());
+ motions.push(xyz.iter().flat_map(|p| [-p[1], p[0]]).collect());
+ }
+ for u in motions {
+ let f = k.mul_vec(&u);
+ let scale = k.data.iter().fold(0.0f64, |m, v| m.max(v.abs()));
+ assert!(f.iter().all(|v| v.abs() < 1e-12 * scale), "{kind:?} {theory:?}: {f:?}");
+ }
+ }
+ }
+ }
+
+ #[test]
+ fn linear_fields_give_exact_constant_stress() {
+ let (young, nu) = (1000.0, 0.3);
+ let d = isotropic(young, nu).unwrap();
+ for (kind, f) in variants() {
+ let xyz = element(kind);
+ // Plane: u = 0.01 x + 0.002 y, v = −0.003 x + 0.004 y.
+ for theory in [Theory::PlaneStress { thickness: 0.7 }, Theory::PlaneStrain { thickness: 1.0 }] {
+ let e = Plane::new(kind, &xyz, f, theory).unwrap();
+ let u: Vec<[f64; 2]> = xyz.iter().map(|p| [0.01 * p[0] + 0.002 * p[1], -0.003 * p[0] + 0.004 * p[1]]).collect();
+ let (exx, eyy, gxy) = (0.01, 0.004, -0.001);
+ let (sxx, syy, szz, ezz) = match theory {
+ Theory::PlaneStress { .. } => {
+ let c = young / (1.0 - nu * nu);
+ (c * (exx + nu * eyy), c * (eyy + nu * exx), 0.0, -nu / (1.0 - nu) * (exx + eyy))
+ }
+ _ => {
+ let l = young * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
+ let m = young / (2.0 * (1.0 + nu));
+ (l * (exx + eyy) + 2.0 * m * exx, l * (exx + eyy) + 2.0 * m * eyy, l * (exx + eyy), 0.0)
+ }
+ };
+ let txy = young / (2.0 * (1.0 + nu)) * gxy;
+ for p in e.stresses(&d, &u).unwrap() {
+ let want = [sxx, syy, szz, txy, 0.0, 0.0];
+ assert!(p.stress.iter().zip(want).all(|(a, b)| (a - b).abs() < 1e-10), "{kind:?} {f:?} {theory:?}: {:?} vs {want:?}", p.stress);
+ assert!((p.strain[2] - ezz).abs() < 1e-14, "{kind:?} {theory:?}");
+ }
+ // The extrapolated nodal values are the constant too.
+ let values: Vec<[f64; 6]> = e.stresses(&d, &u).unwrap().iter().map(|p| p.stress).collect();
+ for v in e.extrapolate(&values).unwrap() {
+ assert!((v[0] - sxx).abs() < 1e-10 && (v[3] - txy).abs() < 1e-10);
+ }
+ }
+ // Axisymmetric: u_r = 0.01 r, u_z = 0.004 z gives εrr = εθθ = 0.01, εzz = 0.004.
+ let e = Plane::new(kind, &xyz, f, Theory::Axisymmetric).unwrap();
+ let u: Vec<[f64; 2]> = xyz.iter().map(|p| [0.01 * p[0], 0.004 * p[1]]).collect();
+ let l = young * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
+ let m = young / (2.0 * (1.0 + nu));
+ let tr = 0.024;
+ let want = [l * tr + 2.0 * m * 0.01, l * tr + 2.0 * m * 0.004, l * tr + 2.0 * m * 0.01, 0.0, 0.0, 0.0];
+ for p in e.stresses(&d, &u).unwrap() {
+ assert!(p.stress.iter().zip(want).all(|(a, b)| (a - b).abs() < 1e-10), "{kind:?} {f:?} axisymmetric: {:?}", p.stress);
+ }
+ // And its stiffness has no out-of-balance force for it inside the element: K u = f
+ // where f are the consistent forces of the boundary tractions; their total along Y
+ // vanishes (no Y body force).
+ let k = e.stiffness(&d).unwrap();
+ let flat: Vec = u.iter().flatten().copied().collect();
+ let fy: f64 = k.mul_vec(&flat).iter().skip(1).step_by(2).sum();
+ assert!(fy.abs() < 1e-10, "{kind:?} {f:?}: {fy}");
+ }
+ }
+
+ #[test]
+ fn volumes_masses_and_loads() {
+ let pi = std::f64::consts::PI;
+ // A unit square at x ∈ [1, 2]: area 1, ring volume 2π · 1.5.
+ let sq = [[1.0, 0.0, 0.0], [2.0, 0.0, 0.0], [2.0, 1.0, 0.0], [1.0, 1.0, 0.0]];
+ let e = Plane::new(ElementKind::Quad4, &sq, Formulation::Full, Theory::PlaneStress { thickness: 3.0 }).unwrap();
+ assert!((e.volume().unwrap() - 3.0).abs() < 1e-14);
+ assert!((e.mass(2.0).unwrap().data.iter().sum::() - 6.0).abs() < 1e-13);
+ let e = Plane::new(ElementKind::Quad4, &sq, Formulation::Full, Theory::Axisymmetric).unwrap();
+ assert!((e.volume().unwrap() - 3.0 * pi).abs() < 1e-13);
+ assert!((e.lumped_mass(1.0).unwrap().iter().sum::() - 3.0 * pi).abs() < 1e-12);
+ // A centrifugal field ω² r: total ∫ ρ ω² r dV = 2π ∫∫ r² dr dz = 2π · 7/3.
+ let f = e.body_force_field(|x| [x[0], 0.0, 0.0]).unwrap();
+ assert!((f.iter().map(|v| v[0]).sum::() - 2.0 * pi * 7.0 / 3.0).abs() < 1e-12);
+ // Pressure on the outer edge (x = 2) of the ring pushes inward with p · 2π · 2 · 1.
+ let fp = modelcraft_model::nodal::edge_pressure(&[sq[1], sq[2]], Theory::Axisymmetric.extent(), 1.0).unwrap();
+ assert!((fp.iter().map(|v| v[0]).sum::() + 4.0 * pi).abs() < 1e-12, "{fp:?}");
+ }
+
+ #[test]
+ fn bad_input_is_an_error() {
+ let d = isotropic(1000.0, 0.3).unwrap();
+ let sq = [[1.0, 0.0, 0.0], [2.0, 0.0, 0.0], [2.0, 1.0, 0.0], [1.0, 1.0, 0.0]];
+ let cw = [sq[0], sq[3], sq[2], sq[1]];
+ let e = Plane::new(ElementKind::Quad4, &cw, Formulation::Full, Theory::PlaneStrain { thickness: 1.0 }).unwrap();
+ assert!(matches!(e.stiffness(&d), Err(Error::Degenerate(_))));
+ assert!(Plane::new(ElementKind::Quad4, &sq, Formulation::Reduced, Theory::Axisymmetric).is_err());
+ assert!(Plane::new(ElementKind::Tri3, &sq[..3], Formulation::IncompatibleModes, Theory::Axisymmetric).is_err());
+ assert!(Plane::new(ElementKind::Hex8, &sq, Formulation::Full, Theory::Axisymmetric).is_err());
+ let off = [[1.0, 0.0, 0.5], sq[1], sq[2], sq[3]];
+ assert!(Plane::new(ElementKind::Quad4, &off, Formulation::Full, Theory::Axisymmetric).is_err());
+ let neg = sq.map(|p| [p[0] - 1.5, p[1], 0.0]);
+ assert!(Plane::new(ElementKind::Quad4, &neg, Formulation::Full, Theory::Axisymmetric).is_err());
+ assert!(Plane::new(ElementKind::Quad4, &sq, Formulation::Full, Theory::PlaneStress { thickness: f64::NAN }).is_err());
+ // A ring element touching the axis is fine (its integration points are off it).
+ let on_axis = sq.map(|p| [p[0] - 1.0, p[1], 0.0]);
+ let e = Plane::new(ElementKind::Quad4, &on_axis, Formulation::IncompatibleModes, Theory::Axisymmetric).unwrap();
+ assert!(e.stiffness(&d).is_ok());
+ assert!(formulation(ElementKind::Quad4, &modelcraft_model::Formulation { hybrid: true, ..Default::default() }).is_err());
+ }
+}
diff --git a/crates/fem/src/quadrature.rs b/crates/fem/src/quadrature.rs
index 44e6790..ff863fa 100644
--- a/crates/fem/src/quadrature.rs
+++ b/crates/fem/src/quadrature.rs
@@ -11,6 +11,7 @@
//! weights only.
//! - On the reference triangle {ξ, η ≥ 0, ξ + η ≤ 1}: the same collapsed product in two
//! dimensions (Jacobian 1 − u), exact for degree 2n − 2.
+//! Also the centroid rule (degree 1) and the symmetric three-point rule (degree 2).
//! - On the reference square [−1, 1]²: tensor Gauss–Legendre.
//! - On the reference wedge (triangle × [−1, 1]): the symmetric three-point triangle rule
//! (points at area coordinates (⅔, ⅙, ⅙) and permutations, degree 2; Zienkiewicz, Taylor &
@@ -57,6 +58,11 @@ pub enum Rule {
Quad(u8),
/// The collapsed product on the reference triangle, n points per direction (z = 0).
TriConical(u8),
+ /// The triangle's centroid, weight 1/2 (degree 1).
+ TriCentroid,
+ /// Three symmetric points on the triangle at area coordinates (⅔, ⅙, ⅙) and permutations,
+ /// weight 1/6 each (degree 2); point k sits near corner k.
+ Tri3,
/// Three triangle points times two Gauss points in ζ on the reference wedge: point k of each
/// layer sits near corner k, the lower layer first.
Wedge,
@@ -109,6 +115,11 @@ impl Rule {
.flat_map(|&(u, wu)| g.iter().map(move |&(v, wv)| Point { xi: [u, v * (1.0 - u), 0.0], weight: wu * wv * (1.0 - u) }))
.collect()
}
+ Rule::TriCentroid => vec![Point { xi: [1.0 / 3.0, 1.0 / 3.0, 0.0], weight: 0.5 }],
+ Rule::Tri3 => {
+ let (a, b) = (2.0 / 3.0, 1.0 / 6.0);
+ [[b, b, 0.0], [a, b, 0.0], [b, a, 0.0]].into_iter().map(|xi| Point { xi, weight: 1.0 / 6.0 }).collect()
+ }
Rule::Wedge => {
let (a, b) = (2.0 / 3.0, 1.0 / 6.0);
let tri = [[b, b], [a, b], [b, a]];
@@ -172,10 +183,12 @@ mod tests {
}
}
// Triangle: ∫ ξ^a η^b = a! b! / (a + b + 2)!.
- for a in 0..=4 {
- for b in 0..=4 - a {
- let got = integrate(Rule::TriConical(3), a, b, 0);
- assert!((got - f(a) * f(b) / f(a + b + 2)).abs() < 1e-14, "tri ξ^{a} η^{b}");
+ for (rule, degree) in [(Rule::TriConical(3), 4), (Rule::TriConical(4), 6), (Rule::TriCentroid, 1), (Rule::Tri3, 2)] {
+ for a in 0..=degree {
+ for b in 0..=degree - a {
+ let got = integrate(rule, a, b, 0);
+ assert!((got - f(a) * f(b) / f(a + b + 2)).abs() < 1e-14, "{rule:?} ξ^{a} η^{b}");
+ }
}
}
}
diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs
index b4e0cd4..7af2384 100644
--- a/crates/mesh/src/lib.rs
+++ b/crates/mesh/src/lib.rs
@@ -9,6 +9,7 @@
pub mod cdt;
pub mod delaunay;
+pub mod planar;
pub mod quad;
pub mod quality;
pub mod size;
diff --git a/crates/mesh/src/planar.rs b/crates/mesh/src/planar.rs
new file mode 100644
index 0000000..0460581
--- /dev/null
+++ b/crates/mesh/src/planar.rs
@@ -0,0 +1,377 @@
+//! 2D models of planar regions: the meshes plane stress, plane strain and axisymmetric analyses
+//! run on.
+//!
+//! A region is bounded by loops of *curves* (polylines; the first loop is the outer boundary, the
+//! others are holes). [`mesh_region`] divides every curve by the size, meshes the region with the
+//! free quadrilateral mesher or the free triangle mesher of [`crate::quad`] on those boundary
+//! nodes, and adds mid-side nodes for quadratic elements (on the boundary polyline for boundary
+//! edges, at the chord midpoint inside). Every element edge on a curve is reported with the
+//! curve's id, so a mesher above can tag it with the geometry it came from.
+//!
+//! [`face_boundary`] gives the curves of one planar face of a tessellated body: the face's
+//! boundary loops in the XY plane, split where the neighbouring face changes (each curve's id is
+//! the neighbour's face index, the B-rep edge the two faces share) and at sharp corners.
+
+use std::collections::BTreeMap;
+
+use modelcraft_model::ElementKind;
+
+use crate::quad::{QuadOptions, free_quad_mesh, free_tri_mesh};
+use crate::size::{UniformSize, divide_polyline, edge_divisions, even_divisions, ideal_divisions};
+use crate::surface::weld::weld;
+use crate::{Error, Result};
+
+/// A boundary curve: a polyline from its start to the start of the next curve of its loop (the
+/// end point is not repeated), and an id reported on the element edges that lie on it.
+#[derive(Debug, Clone, PartialEq)]
+pub struct Curve {
+ pub id: Option,
+ pub points: Vec<[f64; 2]>,
+}
+
+/// A planar mesh of 2D elements in the XY plane.
+#[derive(Debug, Clone, PartialEq)]
+pub struct PlanarMesh {
+ pub kind: ElementKind,
+ pub nodes: Vec<[f64; 2]>,
+ /// Counterclockwise node lists in the model's node order (corners, then mid-side nodes).
+ pub elements: Vec>,
+ /// Element edges on a curve with an id: (element index, edge number from 1, curve id).
+ pub edges: Vec<(usize, u8, u32)>,
+ /// The node at the start of each curve with an id: (node, curve id).
+ pub corners: Vec<(u32, u32)>,
+}
+
+/// Turns sharper than this (radians) start a new curve in [`face_boundary`], so corners stay
+/// nodes of the mesh.
+const CORNER_ANGLE: f64 = std::f64::consts::PI / 6.0;
+
+/// Mesh the region bounded by `loops` with elements of `kind` (tri3, tri6, quad4 or quad8) of
+/// about `size`.
+pub fn mesh_region(loops: &[Vec], kind: ElementKind, size: f64, opts: &QuadOptions) -> Result {
+ if kind.dimension() != 2 {
+ return Err(Error::Invalid(format!("a planar mesh has tri3, tri6, quad4 or quad8 elements, not {}", kind.name())));
+ }
+ if !(size > 0.0 && size.is_finite()) {
+ return Err(Error::Invalid(format!("the element size must be a positive number, got {size}")));
+ }
+ if loops.is_empty() || loops.iter().any(Vec::is_empty) {
+ return Err(Error::Invalid("a region needs at least one loop of curves".into()));
+ }
+ let quads = matches!(kind, ElementKind::Quad4 | ElementKind::Quad8);
+ let quadratic = kind.is_quadratic();
+ let field = UniformSize(size);
+ // Boundary points per loop, the curve of each boundary segment, and the boundary mid points.
+ let mut boundary: Vec> = Vec::with_capacity(loops.len());
+ let mut segment_curve: Vec> = Vec::new();
+ let mut segment_mid: Vec<[f64; 2]> = Vec::new();
+ let mut curve_starts: Vec<(usize, u32)> = Vec::new();
+ for (l, curves) in loops.iter().enumerate() {
+ let polylines: Vec> = curves
+ .iter()
+ .enumerate()
+ .map(|(i, c)| {
+ let next = curves.get((i + 1) % curves.len()).and_then(|n| n.points.first()).copied();
+ c.points.iter().chain(next.as_ref()).map(|p| [p[0], p[1], 0.0]).collect()
+ })
+ .collect();
+ let mut counts: Vec = polylines
+ .iter()
+ .map(|p| if quads { ideal_divisions(&field, p).and_then(even_divisions) } else { edge_divisions(&field, p) })
+ .collect::>()?;
+ // Enough points per loop for the meshers: 3 for triangles; 4 for a lone quad loop, 6 with
+ // holes.
+ let least = if !quads {
+ 3
+ } else if loops.len() > 1 {
+ 6
+ } else {
+ 4
+ };
+ while counts.iter().sum::() < least {
+ let longest = (0..counts.len()).max_by(|&a, &b| counts[a].cmp(&counts[b]).then(b.cmp(&a))).unwrap_or(0);
+ if let Some(c) = counts.get_mut(longest) {
+ *c += if quads { 2 } else { 1 };
+ }
+ }
+ let mut pts = Vec::new();
+ for ((poly, n), curve) in polylines.iter().zip(&counts).zip(curves) {
+ let fine = divide_polyline(&field, poly, if quadratic { 2 * n } else { *n })?;
+ let step = if quadratic { 2 } else { 1 };
+ if let Some(id) = curve.id {
+ curve_starts.push((boundary.iter().map(Vec::len).sum::() + pts.len(), id));
+ }
+ for k in 0..*n {
+ let p = fine.get(step * k).ok_or_else(|| Error::Failed("curve division".into()))?;
+ pts.push([p[0], p[1]]);
+ segment_curve.push(curve.id);
+ let m = if quadratic { fine.get(2 * k + 1) } else { None }.copied().unwrap_or_default();
+ segment_mid.push([m[0], m[1]]);
+ }
+ }
+ if pts.len() < 3 {
+ return Err(Error::Invalid(format!("loop {l} is too short to mesh")));
+ }
+ boundary.push(pts);
+ }
+ let (nodes, corners_of): (Vec<[f64; 2]>, Vec>) = if quads {
+ let m = free_quad_mesh(&boundary, &field, opts)?;
+ (m.nodes, m.quads.iter().map(|q| q.to_vec()).collect())
+ } else {
+ let m = free_tri_mesh(&boundary, &field, opts)?;
+ (m.nodes, m.triangles.iter().map(|t| t.to_vec()).collect())
+ };
+ // Boundary segment of a node pair: segment k of loop l joins its points k and k + 1.
+ let mut segment: BTreeMap<(u32, u32), usize> = BTreeMap::new();
+ let mut offset = 0usize;
+ let mut s = 0usize;
+ for pts in &boundary {
+ let n = pts.len();
+ for k in 0..n {
+ let (a, b) = ((offset + k) as u32, (offset + (k + 1) % n) as u32);
+ segment.insert((a.min(b), a.max(b)), s);
+ s += 1;
+ }
+ offset += n;
+ }
+ let mut nodes = nodes;
+ let mut mids: BTreeMap<(u32, u32), u32> = BTreeMap::new();
+ let mut elements = Vec::with_capacity(corners_of.len());
+ let mut edges = Vec::new();
+ for (e, corners) in corners_of.iter().enumerate() {
+ let mut list = corners.clone();
+ for (i, local) in kind.edges().iter().enumerate() {
+ let (Some(&a), Some(&b)) = (local.first().and_then(|x| corners.get(*x)), local.get(1).and_then(|x| corners.get(*x))) else {
+ return Err(Error::Failed("element edge".into()));
+ };
+ let key = (a.min(b), a.max(b));
+ let on = segment.get(&key).copied();
+ if let Some(id) = on.and_then(|s| segment_curve.get(s).copied().flatten()) {
+ edges.push((e, i as u8 + 1, id));
+ }
+ if quadratic {
+ let id = match mids.get(&key) {
+ Some(id) => *id,
+ None => {
+ let p = match on.and_then(|s| segment_mid.get(s)) {
+ Some(m) => *m,
+ None => {
+ let (pa, pb) =
+ (nodes.get(a as usize).copied().unwrap_or_default(), nodes.get(b as usize).copied().unwrap_or_default());
+ [0.5 * (pa[0] + pb[0]), 0.5 * (pa[1] + pb[1])]
+ }
+ };
+ nodes.push(p);
+ let id = u32::try_from(nodes.len() - 1).map_err(|_| Error::Limit("too many nodes".into()))?;
+ mids.insert(key, id);
+ id
+ }
+ };
+ list.push(id);
+ }
+ }
+ elements.push(list);
+ }
+ let corners = curve_starts.into_iter().map(|(n, id)| (n as u32, id)).collect();
+ Ok(PlanarMesh { kind, nodes, elements, edges, corners })
+}
+
+/// The boundary of one planar face of a triangulated body, in the XY plane.
+#[derive(Debug, Clone, PartialEq)]
+pub struct FaceBoundary {
+ /// The outer loop first, then the holes; each curve's id is the face index of the neighbour
+ /// it borders (`None` where no neighbour was found).
+ pub loops: Vec>,
+ /// The face's Z.
+ pub z: f64,
+}
+
+/// The boundary of face `face` of a triangulated body (`tri_face` gives each triangle's face):
+/// the face must lie in a plane of constant Z (within `tolerance`); points of neighbouring faces
+/// within `tolerance` are the same point.
+pub fn face_boundary(positions: &[[f64; 3]], triangles: &[[u32; 3]], tri_face: &[u32], face: u32, tolerance: f64) -> Result {
+ if !(tolerance >= 0.0 && tolerance.is_finite()) {
+ return Err(Error::Invalid(format!("tolerance must be a non-negative number, got {tolerance}")));
+ }
+ if let Some(v) = triangles.iter().flatten().find(|v| **v as usize >= positions.len()) {
+ return Err(Error::Invalid(format!("triangle vertex {v} is out of range")));
+ }
+ if tri_face.len() != triangles.len() {
+ return Err(Error::Invalid("every triangle needs its face".into()));
+ }
+ // One weld of the whole tessellation: each face's copies of a shared point become one point,
+ // so neighbours meet on welded edges.
+ let welded = weld(positions, triangles, 0.0, tolerance)?;
+ let xyz = welded.points;
+ let face_of = |t: usize| welded.source_triangle.get(t).and_then(|s| tri_face.get(*s as usize)).copied();
+ let mine: Vec<[u32; 3]> = welded.triangles.iter().enumerate().filter(|(t, _)| face_of(*t) == Some(face)).map(|(_, w)| *w).collect();
+ if mine.is_empty() {
+ return Err(Error::Invalid(format!("face {face} has no triangles")));
+ }
+ let z = mine.first().and_then(|w| xyz.get(w[0] as usize)).map_or(0.0, |p| p[2]);
+ let flat = mine.iter().flatten().filter_map(|v| xyz.get(*v as usize)).all(|p| (p[2] - z).abs() <= tolerance.max(1e-12 * (1.0 + z.abs())));
+ if !flat {
+ return Err(Error::Invalid("the face does not lie in a plane parallel to XY (a 2D model lies in the XY plane)".into()));
+ }
+ // Boundary edges: directed edges whose reverse is not in the face.
+ let mut directed: BTreeMap<(u32, u32), usize> = BTreeMap::new();
+ for w in &mine {
+ for k in 0..3 {
+ *directed.entry((w[k], w[(k + 1) % 3])).or_default() += 1;
+ }
+ }
+ let mut next: BTreeMap = BTreeMap::new();
+ for (&(a, b), &count) in &directed {
+ if directed.contains_key(&(b, a)) {
+ continue;
+ }
+ if count > 1 || next.insert(a, b).is_some() {
+ return Err(Error::Invalid("the face's boundary touches itself (a non-manifold face); split the face".into()));
+ }
+ }
+ // Neighbours: the other faces' welded edges, by their ends.
+ let mut other_faces: BTreeMap<(u32, u32), u32> = BTreeMap::new();
+ for (t, w) in welded.triangles.iter().enumerate() {
+ let Some(f) = face_of(t).filter(|f| *f != face) else { continue };
+ for k in 0..3 {
+ let (a, b) = (w[k], w[(k + 1) % 3]);
+ other_faces.entry((a.min(b), a.max(b))).or_insert(f);
+ }
+ }
+ let neighbour = |a: u32, b: u32| other_faces.get(&(a.min(b), a.max(b))).copied();
+ // Chain the boundary into loops, each started at its smallest vertex id.
+ let mut loops: Vec<(f64, Vec)> = Vec::new();
+ let mut unused: BTreeMap = next.clone();
+ while let Some((&start, _)) = unused.iter().next() {
+ let mut chain = vec![start];
+ let mut v = start;
+ loop {
+ let Some(w) = unused.remove(&v) else {
+ return Err(Error::Invalid("the face's boundary is not closed".into()));
+ };
+ if w == start {
+ break;
+ }
+ chain.push(w);
+ v = w;
+ if chain.len() > next.len() {
+ return Err(Error::Failed("boundary chaining did not close".into()));
+ }
+ }
+ let at = |i: usize| chain.get(i % chain.len()).and_then(|v| xyz.get(*v as usize)).copied().unwrap_or_default();
+ let n = chain.len();
+ let ids: Vec> = (0..n).map(|i| neighbour(chain[i], chain[(i + 1) % n])).collect();
+ // A curve starts where the neighbour changes or the boundary turns sharply.
+ let starts: Vec = (0..n)
+ .map(|i| {
+ let (p, q, r) = (at(i + n - 1), at(i), at(i + 1));
+ let (u, w) = ([q[0] - p[0], q[1] - p[1]], [r[0] - q[0], r[1] - q[1]]);
+ let turn = (u[0] * w[1] - u[1] * w[0]).atan2(u[0] * w[0] + u[1] * w[1]).abs();
+ ids[i] != ids[(i + n - 1) % n] || turn > CORNER_ANGLE
+ })
+ .collect();
+ let first = starts.iter().position(|s| *s).unwrap_or(0);
+ let mut curves: Vec = Vec::new();
+ for k in 0..n {
+ let i = (first + k) % n;
+ let p = at(i);
+ if starts[i] || curves.is_empty() {
+ curves.push(Curve { id: ids[i], points: Vec::new() });
+ }
+ if let Some(c) = curves.last_mut() {
+ c.points.push([p[0], p[1]]);
+ }
+ }
+ let area: f64 = (0..n).map(|i| at(i)[0] * at(i + 1)[1] - at(i + 1)[0] * at(i)[1]).sum::() / 2.0;
+ loops.push((area.abs(), curves));
+ }
+ // The outer loop encloses the most area.
+ let outer = (0..loops.len()).max_by(|&a, &b| loops[a].0.total_cmp(&loops[b].0).then(b.cmp(&a))).unwrap_or(0);
+ let mut out = Vec::with_capacity(loops.len());
+ if outer < loops.len() {
+ out.push(loops.remove(outer).1);
+ }
+ out.extend(loops.into_iter().map(|l| l.1));
+ Ok(FaceBoundary { loops: out, z })
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+
+ fn square_with_hole() -> Vec> {
+ let side = |id: u32, a: [f64; 2]| Curve { id: Some(id), points: vec![a] };
+ let outer = vec![side(1, [0.0, 0.0]), side(2, [4.0, 0.0]), side(3, [4.0, 3.0]), side(4, [0.0, 3.0])];
+ let circle: Vec<[f64; 2]> = (0..32).map(|k| 2.0 * std::f64::consts::PI * k as f64 / 32.0).map(|t| [2.0 + t.cos(), 1.5 + t.sin()]).collect();
+ vec![outer, vec![Curve { id: Some(5), points: circle }]]
+ }
+
+ #[test]
+ fn regions_mesh_with_every_element_kind_and_report_curve_edges() {
+ let loops = square_with_hole();
+ let area_want = 12.0 - std::f64::consts::PI;
+ for kind in [ElementKind::Tri3, ElementKind::Tri6, ElementKind::Quad4, ElementKind::Quad8] {
+ let m = mesh_region(&loops, kind, 0.4, &QuadOptions::default()).unwrap();
+ assert!(m.elements.iter().all(|e| e.len() == kind.node_count()));
+ // Area by the corner polygons (the hole's chords make it a little larger).
+ let area: f64 = m
+ .elements
+ .iter()
+ .map(|e| {
+ let c: Vec<[f64; 2]> = e[..kind.corner_count()].iter().map(|v| m.nodes[*v as usize]).collect();
+ (0..c.len()).map(|i| c[i][0] * c[(i + 1) % c.len()][1] - c[(i + 1) % c.len()][0] * c[i][1]).sum::() / 2.0
+ })
+ .sum();
+ assert!(area > area_want && area - area_want < 0.1, "{kind:?}: {area}");
+ // Edges on each curve: the right side (x = 4, curve 2) has length 3.
+ let on_right: Vec<_> = m.edges.iter().filter(|e| e.2 == 2).collect();
+ assert!(!on_right.is_empty());
+ for (e, edge, _) in on_right {
+ let local = kind.edges()[usize::from(*edge) - 1];
+ for i in local {
+ assert!((m.nodes[m.elements[*e][*i] as usize][0] - 4.0).abs() < 1e-12, "{kind:?}");
+ }
+ }
+ // Quadratic mid-side nodes on the hole lie on its polyline (radius 1 within the chord sag).
+ assert_eq!(m.corners.len(), 5);
+ assert_eq!(mesh_region(&loops, kind, 0.4, &QuadOptions::default()).unwrap(), m, "deterministic");
+ }
+ assert!(mesh_region(&loops, ElementKind::Hex8, 0.4, &QuadOptions::default()).is_err());
+ assert!(mesh_region(&loops, ElementKind::Quad4, 0.0, &QuadOptions::default()).is_err());
+ assert!(mesh_region(&loops, ElementKind::Quad4, f64::NAN, &QuadOptions::default()).is_err());
+ assert!(mesh_region(&[], ElementKind::Quad4, 1.0, &QuadOptions::default()).is_err());
+ assert!(mesh_region(&loops, ElementKind::Tri3, 0.02, &QuadOptions { max_elements: 1000, ..QuadOptions::default() }).is_err());
+ }
+
+ #[test]
+ fn a_tessellated_face_gives_its_loops_and_neighbours() {
+ // A unit cube's bottom face (z = 0, faces 0), split into two triangles; four side faces 1–4
+ // with their own copies of the points.
+ let p = vec![
+ [0.0, 0.0, 0.0],
+ [1.0, 0.0, 0.0],
+ [1.0, 1.0, 0.0],
+ [0.0, 1.0, 0.0],
+ // Side y = 0 (face 1).
+ [0.0, 0.0, 0.0],
+ [1.0, 0.0, 0.0],
+ [1.0, 0.0, 1.0],
+ // Side x = 1 (face 2).
+ [1.0, 0.0, 0.0],
+ [1.0, 1.0, 0.0],
+ [1.0, 1.0, 1.0],
+ ];
+ let tris = vec![[0, 2, 1], [0, 3, 2], [4, 5, 6], [7, 8, 9]];
+ let faces = vec![0, 0, 1, 2];
+ let b = face_boundary(&p, &tris, &faces, 0, 1e-9).unwrap();
+ assert_eq!(b.loops.len(), 1);
+ let ids: Vec> = b.loops[0].iter().map(|c| c.id).collect();
+ assert!(ids.contains(&Some(1)) && ids.contains(&Some(2)) && ids.contains(&None), "{ids:?}");
+ let points: usize = b.loops[0].iter().map(|c| c.points.len()).sum();
+ assert_eq!(points, 4);
+ assert!(face_boundary(&p, &tris, &faces, 1, 1e-9).is_err(), "face 1 is vertical");
+ assert!(face_boundary(&p, &tris, &faces, 7, 1e-9).is_err());
+ assert!(face_boundary(&p, &[[0, 1, 99]], &[0], 0, 1e-9).is_err());
+ assert!(face_boundary(&p, &tris, &faces, 0, f64::NAN).is_err());
+ }
+}
diff --git a/crates/mesh/src/quad/mod.rs b/crates/mesh/src/quad/mod.rs
index 313787e..2a52362 100644
--- a/crates/mesh/src/quad/mod.rs
+++ b/crates/mesh/src/quad/mod.rs
@@ -160,6 +160,51 @@ pub fn free_quad_mesh(loops: &[Vec<[f64; 2]>], size: &dyn SizeField, opts: &Quad
free::mesh(loops, size, opts)
}
+/// A planar triangle mesh.
+#[derive(Debug, Clone, PartialEq)]
+pub struct TriMesh {
+ /// Node coordinates; `nodes[..boundary_nodes]` are the input boundary points in input order.
+ pub nodes: Vec<[f64; 2]>,
+ /// Counterclockwise corner node ids (the `Tri3` node order).
+ pub triangles: Vec<[u32; 3]>,
+ pub boundary_nodes: usize,
+}
+
+/// Free triangle mesh of a planar region on prescribed boundary nodes: the constrained Delaunay
+/// triangulation of the loops ([`cdt`]: Bowyer–Watson, Sloan 1993 recovery), refined to the size
+/// field by Ruppert-style circumcentre insertion that never splits a boundary edge (Steiner points
+/// keep a clearance of half the local boundary spacing from the boundary), so the boundary nodes
+/// are exactly the input points. Every triangle is checked counterclockwise with exact
+/// predicates. Loops as for [`free_quad_mesh`], with any number of points (at least 3 each).
+pub fn free_tri_mesh(loops: &[Vec<[f64; 2]>], size: &dyn SizeField, opts: &QuadOptions) -> Result {
+ let region = domain::region(loops)?;
+ let mut chords = Vec::new();
+ for l in ®ion.loops {
+ for (i, &a) in l.iter().enumerate() {
+ let b = l.get((i + 1) % l.len()).copied().ok_or_else(|| Error::Failed("loop index".into()))?;
+ let (pa, pb) = (*region.point(a)?, *region.point(b)?);
+ let (ca, cb) = (pa.coords(), pb.coords());
+ let m = Point2::new(snap_to_domain(0.5 * (ca[0] + cb[0])), snap_to_domain(0.5 * (ca[1] + cb[1])))?;
+ chords.push((pa, pb, m));
+ }
+ }
+ let guard = guard::Guard::new(&chords)?;
+ let mut cdt = cdt::Cdt::build(®ion, opts.max_elements.saturating_add(4 * region.pts.len() + 16))?;
+ cdt.refine(size, 1.0, &|p, h| guard.allows(p, cdt::CLEARANCE * h))?;
+ let tri = cdt.into_triangulation();
+ if tri.tris.len() > opts.max_elements {
+ return Err(Error::Limit(format!("the mesh would have {} triangles; the limit is {}", tri.tris.len(), opts.max_elements)));
+ }
+ for t in &tri.tris {
+ let [a, b, c] = t.map(|v| tri.p2.get(v as usize));
+ match (a, b, c) {
+ (Some(a), Some(b), Some(c)) if crate::predicates::orient2d(a, b, c) == Sign::Positive => {}
+ _ => return Err(Error::Failed("a triangle is inverted or degenerate".into())),
+ }
+ }
+ Ok(TriMesh { nodes: tri.pts, triangles: tri.tris, boundary_nodes: region.pts.len() })
+}
+
/// Mapped mesh of a four-sided region: `boundary` lists the region's boundary nodes once (either
/// orientation), and `corners` are the positions in it of the four corners in boundary order.
/// Opposite sides need the same number of divisions. The grid is transfinite interpolation
diff --git a/crates/mesh/src/quad/tests.rs b/crates/mesh/src/quad/tests.rs
index e9b7942..ffb0e12 100644
--- a/crates/mesh/src/quad/tests.rs
+++ b/crates/mesh/src/quad/tests.rs
@@ -537,3 +537,34 @@ fn unsmoothed_meshes_are_valid() {
let m = free_quad_mesh(&loops, &UniformSize(0.2), &QuadOptions { smoothing_sweeps: 0, ..opts() }).unwrap();
check(&m, &loops);
}
+
+#[test]
+fn free_triangle_mesh_of_a_plate_with_a_hole() {
+ // Odd point counts are fine for triangles.
+ let outer = rect(0.0, 0.0, 10.0, 6.0, 9, 5);
+ let hole: Vec<[f64; 2]> = circle(4.0, 3.0, 1.5, 17).into_iter().rev().collect();
+ let loops = vec![outer.clone(), hole.clone()];
+ let m = free_tri_mesh(&loops, &UniformSize(0.7), &opts()).unwrap();
+ assert_eq!(m.boundary_nodes, outer.len() + hole.len());
+ let boundary: Vec<[f64; 2]> = outer.iter().chain(&hole).copied().collect();
+ assert_eq!(&m.nodes[..m.boundary_nodes], boundary.as_slice(), "boundary nodes are the input points");
+ let area: f64 = m
+ .triangles
+ .iter()
+ .map(|t| {
+ let [a, b, c] = t.map(|v| m.nodes[v as usize]);
+ let s = 0.5 * ((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]));
+ assert!(s > 0.0, "counterclockwise");
+ s
+ })
+ .sum();
+ let want = signed_area(&outer).abs() - signed_area(&hole).abs();
+ assert!((area - want).abs() < 1e-9 * want, "{area} vs {want}");
+ assert!(m.triangles.len() > 100, "refined to the size: {}", m.triangles.len());
+ // Deterministic.
+ assert_eq!(free_tri_mesh(&loops, &UniformSize(0.7), &opts()).unwrap(), m);
+ // Hostile.
+ assert!(free_tri_mesh(&[vec![[0.0, 0.0], [1.0, 0.0]]], &UniformSize(0.5), &opts()).is_err());
+ assert!(free_tri_mesh(std::slice::from_ref(&outer), &UniformSize(0.001), &QuadOptions { max_elements: 100, ..opts() }).is_err());
+ assert!(free_tri_mesh(&[vec![[0.0, 0.0], [f64::NAN, 0.0], [0.0, 1.0]]], &UniformSize(0.5), &opts()).is_err());
+}
diff --git a/crates/mesh/src/surface/mod.rs b/crates/mesh/src/surface/mod.rs
index 462c911..6e6d3df 100644
--- a/crates/mesh/src/surface/mod.rs
+++ b/crates/mesh/src/surface/mod.rs
@@ -46,7 +46,7 @@
mod bvh;
pub mod intersect;
mod topology;
-mod weld;
+pub(crate) mod weld;
#[cfg(test)]
mod tests;
diff --git a/crates/model/src/check.rs b/crates/model/src/check.rs
index 651a343..4146575 100644
--- a/crates/model/src/check.rs
+++ b/crates/model/src/check.rs
@@ -10,7 +10,7 @@ use std::collections::{BTreeMap, BTreeSet};
use serde::{Deserialize, Serialize};
use crate::loads::Distribution;
-use crate::{BoundaryCondition, Elastic, Error, Load, Model, Procedure, Result, SectionKind, Target};
+use crate::{BoundaryCondition, Elastic, ElementKind, Error, Extent, Load, Model, Procedure, Result, SectionKind, Target};
/// How serious a [`Diagnostic`] is.
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)]
@@ -47,11 +47,14 @@ pub struct ResolvedSet {
pub nodes: Vec,
pub elements: Vec,
pub faces: Vec,
+ /// Element edges (every element kind; 2D models load the edges of their 2D elements).
+ #[serde(default)]
+ pub edges: Vec,
}
impl ResolvedSet {
pub fn is_empty(&self) -> bool {
- self.nodes.is_empty() && self.elements.is_empty() && self.faces.is_empty()
+ self.nodes.is_empty() && self.elements.is_empty() && self.faces.is_empty() && self.edges.is_empty()
}
}
@@ -70,9 +73,30 @@ pub struct Resolved {
pub nodal_loads: BTreeMap>,
/// The check's warnings (a resolved model has no errors).
pub warnings: Vec,
+ /// How each element of a 2D section (plane stress, plane strain, axisymmetric) extends out
+ /// of its plane; empty for a model without 2D sections. A model is 2D or not as a whole.
+ pub planar: BTreeMap,
}
impl Resolved {
+ /// Whether the model is a 2D model (its sections are plane stress, plane strain or
+ /// axisymmetric).
+ pub fn is_planar(&self) -> bool {
+ !self.planar.is_empty()
+ }
+
+ /// How a 2D model extends out of its plane: revolved when any of its sections is
+ /// axisymmetric (the check refuses a mix), else the first element's thickness; `None` for a
+ /// model without 2D sections.
+ pub fn extent(&self) -> Option {
+ if self.planar.values().any(|e| *e == Extent::Revolved) { Some(Extent::Revolved) } else { self.planar.values().next().copied() }
+ }
+
+ /// The edges of a set that belong to 2D elements: where a 2D model's surface loads act.
+ pub fn planar_edges(&self, set: &str) -> Vec {
+ self.sets.get(set).map(|s| s.edges.iter().filter(|e| self.planar.contains_key(&e.element)).copied().collect()).unwrap_or_default()
+ }
+
/// The nodes a boundary condition, force or moment acts on.
pub fn target_nodes(&self, target: &Target) -> &[u32] {
match target {
@@ -160,6 +184,14 @@ impl Model {
}
resolved.sets.insert(name.clone(), r);
}
+ // The elements of 2D sections.
+ for sec in self.sections.values() {
+ if let (Some(extent), Some(r)) = (sec.kind.extent(), resolved.sets.get(&sec.set)) {
+ for &e in &r.elements {
+ resolved.planar.insert(e, extent);
+ }
+ }
+ }
let empty = crate::ResolvedSet::default();
let set_of = |name: &str| resolved.sets.get(name).unwrap_or(&empty);
@@ -191,12 +223,16 @@ impl Model {
if r.elements.is_empty() {
d.error(&subject, format!("set `{}` has no elements", sec.set));
}
- let fits = |k: crate::ElementKind| match sec.kind {
+ let fits = |k: 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,
+ SectionKind::Shell { .. } | SectionKind::PlaneStress { .. } | SectionKind::PlaneStrain { .. } | SectionKind::Axisymmetric {} => {
+ k.dimension() == 2
+ }
+ SectionKind::Beam { .. } => k == ElementKind::Beam2,
+ SectionKind::Rod { .. } => k == ElementKind::Rod2,
};
+ // The linear element with an incompatible-modes formulation for this kind of section.
+ let incompatible_kind = if sec.kind.is_planar() { ElementKind::Quad4 } else { ElementKind::Hex8 };
let mut wrong = Vec::new();
let mut not_hex8 = Vec::new();
let mut along = Vec::new();
@@ -206,7 +242,7 @@ impl Model {
if !fits(el.kind) {
wrong.push(e);
}
- if sec.formulation.incompatible_modes && el.kind != crate::ElementKind::Hex8 {
+ if sec.formulation.incompatible_modes && el.kind != incompatible_kind {
not_hex8.push(e);
}
if let SectionKind::Beam { orientation, .. } = sec.kind
@@ -225,8 +261,17 @@ impl Model {
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)));
+ d.error(
+ &subject,
+ format!(
+ "incompatible modes apply to {} elements only, not elements {} (in a {} section)",
+ incompatible_kind.name(),
+ first_labels(¬_hex8),
+ sec.kind.name()
+ ),
+ );
}
+
if !along.is_empty() {
d.error(&subject, format!("the orientation vector is parallel to elements {}", first_labels(&along)));
}
@@ -241,6 +286,7 @@ impl Model {
if !unsectioned.is_empty() {
d.error("model", format!("elements {} have no section", first_labels(&unsectioned)));
}
+ self.check_planar(&resolved, &mut d);
// Mass: frequency steps and gravity need a density in every material in use.
let needs_mass = self
@@ -254,7 +300,13 @@ impl Model {
};
Some(format!("step `{}` {why}", s.name))
})
- .or_else(|| self.loads.iter().find(|(_, l)| matches!(l, Load::Gravity { .. })).map(|(n, _)| format!("load `{n}` is gravity")));
+ .or_else(|| {
+ self.loads.iter().find_map(|(n, l)| match l {
+ Load::Gravity { .. } => Some(format!("load `{n}` is gravity")),
+ Load::Centrifugal { .. } => Some(format!("load `{n}` is centrifugal")),
+ _ => None,
+ })
+ });
if let Some(why) = needs_mass {
let used: BTreeSet<&str> = self.sections.values().map(|s| s.material.as_str()).collect();
for m in used {
@@ -285,9 +337,9 @@ impl Model {
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())
+ .iter()
+ .filter(|(l, e)| (e.kind.dimension() == 2 && !resolved.planar.contains_key(l)) || e.kind == ElementKind::Beam2)
+ .flat_map(|(_, e)| e.nodes.iter().copied())
.collect()
});
set_of(s).nodes.iter().any(|n| rotational.contains(n))
@@ -325,8 +377,8 @@ impl Model {
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::Pressure { set, .. } if set_of(set).faces.is_empty() && resolved.planar_edges(set).is_empty() => {
+ d.error(&subject, format!("set `{set}` has no element faces (or 2D element edges) 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 } => {
@@ -338,8 +390,12 @@ impl Model {
d.error(&subject, format!("set `{s}` has no elements"));
}
}
+ Load::Centrifugal { set: Some(s), .. } if set_of(s).elements.is_empty() => d.error(&subject, format!("set `{s}` has no elements")),
_ => {}
}
+ if resolved.is_planar() {
+ Self::check_planar_load(load, resolved.extent(), &subject, &mut d);
+ }
}
// Steps and the activity table.
@@ -406,11 +462,12 @@ impl Model {
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 edges = r.planar_edges(set);
+ if faces.is_empty() && edges.is_empty() {
+ return Err(format!("set `{set}` has no element faces (or 2D element edges) 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()));
+ let w = if faces.is_empty() { self.edge_weights(&edges, &r.planar) } else { self.surface_weights(faces) };
+ return Ok(Some(w.map_err(|e| e.to_string())?.into_iter().map(|(n, a)| (n, scale(*vector, a))).collect()));
}
_ => return Ok(None),
};
@@ -429,14 +486,16 @@ impl Model {
if set.nodes.is_empty() {
return Err(format!("set `{name}` has no nodes"));
}
+ let edges = if set.faces.is_empty() { r.planar_edges(name) } else { Vec::new() };
Ok(Some(match distribution {
Distribution::PerNode => set.nodes.iter().map(|n| (*n, vector)).collect(),
- Distribution::Total if set.faces.is_empty() => {
+ Distribution::Total if set.faces.is_empty() && edges.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 w = if edges.is_empty() { self.surface_weights(&set.faces) } else { self.edge_weights(&edges, &r.planar) };
+ let w = w.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"));
@@ -446,6 +505,92 @@ impl Model {
}))
}
+ /// The rules of a 2D model: one kind of 2D model (plane or axisymmetric) and nothing else, in
+ /// the XY plane (x ≥ 0 when axisymmetric), elements numbered counter-clockwise.
+ fn check_planar(&self, r: &Resolved, d: &mut Diagnostics) {
+ if !r.is_planar() {
+ return;
+ }
+ let others: Vec<&str> = self.sections.iter().filter(|(_, s)| !s.kind.is_planar()).map(|(n, _)| n.as_str()).collect();
+ if !others.is_empty() {
+ d.error(
+ "model",
+ format!(
+ "2D sections (plane stress, plane strain, axisymmetric) and other sections ({}) cannot be mixed in one model yet",
+ others.join(", ")
+ ),
+ );
+ }
+ let axisymmetric = r.extent() == Some(Extent::Revolved);
+ if axisymmetric && r.planar.values().any(|e| *e != Extent::Revolved) {
+ d.error("model", "axisymmetric and plane (stress or strain) sections cannot be mixed in one model");
+ }
+ let (mut off_plane, mut negative_radius, mut clockwise) = (Vec::new(), Vec::new(), Vec::new());
+ for &label in r.planar.keys() {
+ let Some(e) = self.elements.get(&label) else { continue };
+ let xyz: Vec<[f64; 3]> = e.nodes.iter().filter_map(|n| self.nodes.get(n)).map(|n| n.xyz).collect();
+ if xyz.iter().any(|p| p[2] != 0.0) {
+ off_plane.push(label);
+ }
+ if axisymmetric && xyz.iter().any(|p| p[0] < 0.0) {
+ negative_radius.push(label);
+ }
+ // Twice the signed area of the corner polygon (shoelace).
+ let corners = xyz.get(..e.kind.corner_count()).unwrap_or_default();
+ let area2: f64 = corners.iter().zip(corners.iter().cycle().skip(1)).map(|(a, b)| a[0] * b[1] - b[0] * a[1]).sum();
+ if !(area2 > 0.0) {
+ clockwise.push(label);
+ }
+ }
+ if !off_plane.is_empty() {
+ d.error("model", format!("2D elements {} have nodes off the XY plane (z ≠ 0); 2D models lie in the XY plane", first_labels(&off_plane)));
+ }
+ if !negative_radius.is_empty() {
+ d.error(
+ "model",
+ format!(
+ "axisymmetric elements {} have nodes at negative x; x is the radius, so the model lies in x ≥ 0",
+ first_labels(&negative_radius)
+ ),
+ );
+ }
+ if !clockwise.is_empty() {
+ d.error(
+ "model",
+ format!(
+ "2D elements {} are clockwise or have no area in the XY plane (2D elements number their nodes counter-clockwise)",
+ first_labels(&clockwise)
+ ),
+ );
+ }
+ }
+
+ /// A load in a 2D model: in the XY plane, and a rotation about the axis the model can carry.
+ fn check_planar_load(load: &Load, extent: Option, subject: &str, d: &mut Diagnostics) {
+ let in_plane = match load {
+ Load::Force { vector, .. } | Load::Traction { vector, .. } | Load::BodyForce { vector, .. } => vector[2] == 0.0,
+ Load::Gravity { acceleration, .. } => acceleration[2] == 0.0,
+ Load::Moment { .. } | Load::Pressure { .. } => true,
+ Load::Centrifugal { point, axis, .. } => {
+ let ok = if extent == Some(Extent::Revolved) {
+ axis[0] == 0.0 && axis[2] == 0.0 && point[0] == 0.0 && point[2] == 0.0
+ } else {
+ axis[0] == 0.0 && axis[1] == 0.0
+ };
+ if !ok {
+ d.error(
+ subject,
+ "in a 2D model a rotation is about the Y axis (axisymmetric: through x = 0, z = 0) or about an axis along Z (plane models)",
+ );
+ }
+ true
+ }
+ };
+ if !in_plane {
+ d.error(subject, "a 2D model is loaded in the XY plane; this load has a Z component");
+ }
+ }
+
/// 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 {
diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs
index 652a737..9d19fc0 100644
--- a/crates/model/src/lib.rs
+++ b/crates/model/src/lib.rs
@@ -53,10 +53,10 @@ pub use check::{Diagnostic, Resolved, ResolvedSet, Severity};
pub use constraints::{Coupling, CouplingKind, ReferencePoint};
pub use fluids::{BoxFace, BoxFaces, FaceCondition, FluidBoundary, FluidDomain, Profile, SurfaceBoundary, SurfaceCondition};
pub use geometry::{Classification, Entity, GeometryRef};
-pub use loads::{Distribution, Load};
+pub use loads::{Distribution, Load, centrifugal_acceleration};
pub use materials::{Elastic, Fluid, Material};
-pub use mesh::{Element, ElementFace, ElementKind, Node, check_coordinates};
-pub use sections::{BeamProfile, Formulation, Integration, Section, SectionKind};
+pub use mesh::{Element, ElementEdge, ElementFace, ElementKind, Node, check_coordinates};
+pub use sections::{BeamProfile, Extent, Formulation, Integration, Section, SectionKind};
pub use sets::Set;
pub use steps::{BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, InitialVelocity, OutputVariable, Procedure, RigidWall, Step};
pub use units::UnitSystem;
diff --git a/crates/model/src/loads.rs b/crates/model/src/loads.rs
index 4846d55..d14da26 100644
--- a/crates/model/src/loads.rs
+++ b/crates/model/src/loads.rs
@@ -1,4 +1,9 @@
//! Loads. Which steps a load is active in is listed explicitly by each step ([`crate::steps`]).
+//!
+//! In a 2D model (plane stress, plane strain, axisymmetric sections) loads lie in the XY plane:
+//! surface loads act on the element edges of the set, forces are totals over the thickness (or
+//! over the full circumference for axisymmetric models), and a vector with a Z component is an
+//! error.
use std::collections::BTreeMap;
@@ -55,6 +60,27 @@ pub enum Load {
},
/// A force per volume on the set's elements.
BodyForce { set: String, vector: [f64; 3] },
+ /// The centrifugal load of a steady rotation at `angular_velocity` (radians per time unit)
+ /// about the axis through `point` along `axis`: a force per volume ρ ω² r on the mass of the
+ /// set's elements (or the whole model), r the distance from the axis, pointing away from it.
+ Centrifugal {
+ angular_velocity: f64,
+ point: [f64; 3],
+ axis: [f64; 3],
+ #[serde(default, skip_serializing_if = "Option::is_none")]
+ set: Option,
+ },
+}
+
+/// The centrifugal acceleration ω² r at `x` of a rotation at `angular_velocity` about the axis
+/// through `point` along `axis` (any length, not zero): r is the part of `x − point` normal to
+/// the axis, so the acceleration points away from the axis.
+pub fn centrifugal_acceleration(angular_velocity: f64, point: [f64; 3], axis: [f64; 3], x: [f64; 3]) -> [f64; 3] {
+ let len2 = axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2];
+ let d = [x[0] - point[0], x[1] - point[1], x[2] - point[2]];
+ let along = if len2 > 0.0 { (d[0] * axis[0] + d[1] * axis[1] + d[2] * axis[2]) / len2 } else { 0.0 };
+ let w2 = angular_velocity * angular_velocity;
+ [0, 1, 2].map(|k| w2 * (d[k] - along * axis[k]))
}
impl Load {
@@ -67,6 +93,7 @@ impl Load {
Load::Pressure { .. } => "pressure",
Load::Gravity { .. } => "gravity",
Load::BodyForce { .. } => "body_force",
+ Load::Centrifugal { .. } => "centrifugal",
}
}
@@ -75,7 +102,7 @@ impl Load {
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(),
+ Load::Gravity { set, .. } | Load::Centrifugal { set, .. } => set.as_deref(),
}
}
@@ -95,6 +122,21 @@ impl Load {
finite("pressure", *magnitude)?;
return Ok(());
}
+ Load::Centrifugal { angular_velocity, point, axis, .. } => {
+ let w = finite("angular velocity", *angular_velocity)?;
+ if w == 0.0 || !(w * w).is_finite() {
+ return Err(Error::Invalid(format!(
+ "a centrifugal load's angular velocity must be a non-zero number whose square is finite, got {w}"
+ )));
+ }
+ for c in point.iter().chain(axis) {
+ finite("centrifugal axis", *c)?;
+ }
+ if axis.iter().all(|c| *c == 0.0) {
+ return Err(Error::Invalid("a centrifugal load's axis direction must not be zero".into()));
+ }
+ return Ok(());
+ }
};
for c in vector {
finite(self.kind(), *c)?;
diff --git a/crates/model/src/mesh.rs b/crates/model/src/mesh.rs
index 75428c4..019ae13 100644
--- a/crates/model/src/mesh.rs
+++ b/crates/model/src/mesh.rs
@@ -195,6 +195,16 @@ pub struct ElementFace {
pub face: u8,
}
+/// One edge of an element: the element's label and its edge number (1-based,
+/// [`ElementKind::edges`]). In a 2D model (plane stress, plane strain, axisymmetric) the edges on
+/// the boundary are what pressures and tractions act on, as faces are for solids.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)]
+#[serde(deny_unknown_fields)]
+pub struct ElementEdge {
+ pub element: u32,
+ pub edge: u8,
+}
+
/// An element: a topology and its node labels.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct Element {
diff --git a/crates/model/src/nodal.rs b/crates/model/src/nodal.rs
index 1a22dc0..399b50b 100644
--- a/crates/model/src/nodal.rs
+++ b/crates/model/src/nodal.rs
@@ -10,7 +10,7 @@
use std::collections::BTreeMap;
-use crate::{ElementFace, Error, Model, Result};
+use crate::{ElementEdge, ElementFace, Error, Extent, Model, Result};
/// Shape functions and their derivatives at one point of a reference face.
struct Shape {
@@ -137,7 +137,82 @@ pub fn face_normal_weights(xyz: &[[f64; 3]]) -> Result> {
Ok(out)
}
+/// Visit the 3-point Gauss–Legendre points of an element edge (its ends, then the mid-side node
+/// of a quadratic edge): the shape function values, the position and the weighted tangent
+/// (∂x/∂s) w, whose length is the point's share of the edge's length.
+fn integrate_edge(xyz: &[[f64; 3]], mut visit: impl FnMut(&[f64], [f64; 3], [f64; 3])) -> Result<()> {
+ if xyz.iter().flatten().any(|c| !c.is_finite()) {
+ return Err(Error::Invalid("edge coordinates must be finite".into()));
+ }
+ if !matches!(xyz.len(), 2 | 3) {
+ return Err(Error::Invalid(format!("an edge has 2 or 3 nodes, not {}", xyz.len())));
+ }
+ let g = (0.6f64).sqrt();
+ for (s, w) in [(-g, 5.0 / 9.0), (0.0, 8.0 / 9.0), (g, 5.0 / 9.0)] {
+ let (n, dn) = if xyz.len() == 2 {
+ ([(1.0 - s) / 2.0, (1.0 + s) / 2.0, 0.0], [-0.5, 0.5, 0.0])
+ } else {
+ ([s * (s - 1.0) / 2.0, s * (s + 1.0) / 2.0, 1.0 - s * s], [s - 0.5, s + 0.5, -2.0 * s])
+ };
+ let (mut x, mut t) = ([0.0; 3], [0.0; 3]);
+ for ((p, ni), di) in xyz.iter().zip(&n).zip(&dn) {
+ for k in 0..3 {
+ x[k] += ni * p[k];
+ t[k] += di * p[k];
+ }
+ }
+ visit(n.get(..xyz.len()).unwrap_or_default(), x, t.map(|v| v * w));
+ }
+ Ok(())
+}
+
+/// ∫ Nᵢ n e dL for every node of an element edge in the XY plane (node coordinates in the edge's
+/// order), with n the normal to the right of the edge's direction (outward for an edge of a
+/// counter-clockwise 2D element) and e the out-of-plane extent (thickness, or 2πr when
+/// revolved): the consistent nodal forces of a unit pressure pulling outward.
+pub fn edge_normal_weights(xyz: &[[f64; 3]], extent: Extent) -> Result> {
+ let mut out = vec![[0.0; 3]; xyz.len()];
+ integrate_edge(xyz, |n, x, t| {
+ let e = extent.at(x);
+ for (o, ni) in out.iter_mut().zip(n) {
+ o[0] += ni * t[1] * e;
+ o[1] -= ni * t[0] * e;
+ }
+ })?;
+ Ok(out)
+}
+
+/// Consistent nodal forces of a uniform pressure on an element edge in the XY plane (positive
+/// pushes into a counter-clockwise element): `−p` times [`edge_normal_weights`].
+pub fn edge_pressure(xyz: &[[f64; 3]], extent: Extent, magnitude: f64) -> Result> {
+ Ok(edge_normal_weights(xyz, extent)?.into_iter().map(|w| w.map(|c| -magnitude * c)).collect())
+}
+
impl Model {
+ /// The nodal weights `∫ N_i e dL` of a set of 2D element edges (e the out-of-plane extent of
+ /// each element in `planar`), summed per node. Their sum is the area of the surface the
+ /// edges stand for: length times thickness, or the surface of revolution.
+ pub fn edge_weights(&self, edges: &[ElementEdge], planar: &BTreeMap) -> Result> {
+ let mut out: BTreeMap = BTreeMap::new();
+ for &g in edges {
+ let extent = planar.get(&g.element).copied().ok_or_else(|| Error::Invalid(format!("element {} is not a 2D element", g.element)))?;
+ let labels = self.edge_nodes(g).ok_or_else(|| Error::NotFound(format!("element {} edge {}", g.element, g.edge)))?;
+ 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::>()?;
+ integrate_edge(&xyz, |n, x, t| {
+ let dl = (t[0] * t[0] + t[1] * t[1] + t[2] * t[2]).sqrt() * extent.at(x);
+ for (label, ni) in labels.iter().zip(n) {
+ *out.entry(*label).or_default() += ni * dl;
+ }
+ })
+ .map_err(|e| match e {
+ Error::Invalid(m) => Error::Invalid(format!("element {} edge {}: {m}", g.element, g.edge)),
+ e => e,
+ })?;
+ }
+ Ok(out)
+ }
+
/// 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> {
@@ -165,6 +240,7 @@ impl Model {
mod tests {
use super::*;
use crate::ElementKind;
+ use std::collections::BTreeMap;
/// A flat 2 × 3 face of a hex8 or hex20 (face 1, z = 0), and its weights.
fn face_weights(kind: ElementKind) -> Vec {
@@ -211,4 +287,30 @@ mod tests {
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:?}");
}
+
+ #[test]
+ fn edge_weights_are_the_consistent_nodal_loads() {
+ let close = |a: &[f64], b: &[f64]| a.iter().zip(b).all(|(x, y)| (x - y).abs() < 1e-12);
+ // A quad8 2 × 1 in the XY plane: edge 1 (y = 0, x from 0 to 2) of thickness 3 shares
+ // its area 6 as 1/6, 1/6, 2/3; revolved about Y it shares the annulus π(2² − 0) = 4π
+ // as ∫ N 2π x dx.
+ let mut m = Model::default();
+ let p = [[0., 0., 0.], [2., 0., 0.], [2., 1., 0.], [0., 1., 0.], [1., 0., 0.], [2., 0.5, 0.], [1., 1., 0.], [0., 0.5, 0.]];
+ for (i, q) in p.iter().enumerate() {
+ m.set_node(i as u32 + 1, *q).unwrap();
+ }
+ m.set_element(1, ElementKind::Quad8, (1..=8).collect()).unwrap();
+ let edge = [ElementEdge { element: 1, edge: 1 }];
+ let w = m.edge_weights(&edge, &BTreeMap::from([(1, Extent::Thickness(3.0))])).unwrap();
+ assert!(close(&[w[&1], w[&2], w[&5]], &[1.0, 1.0, 4.0]), "{w:?}");
+ let w = m.edge_weights(&edge, &BTreeMap::from([(1, Extent::Revolved)])).unwrap();
+ // ∫₀² N 2π x dx with N on [0, 2]: corner at 0 → 0, corner at 2 → 2π·2/3·... summed: 4π.
+ let pi = std::f64::consts::PI;
+ assert!(close(&[w[&1], w[&2], w[&5]], &[0.0, 4.0 * pi / 3.0, 8.0 * pi / 3.0]), "{w:?}");
+ assert!(m.edge_weights(&edge, &BTreeMap::new()).is_err(), "not a 2D element");
+ // The outward normal of a counter-clockwise element's bottom edge is −Y.
+ let n = edge_normal_weights(&[p[0], p[1], p[4]], Extent::Thickness(1.0)).unwrap();
+ assert!(close(&n.iter().map(|v| v[1]).collect::>(), &[-1.0 / 3.0, -1.0 / 3.0, -4.0 / 3.0]) && n.iter().all(|v| v[0] == 0.0), "{n:?}");
+ assert!(edge_normal_weights(&[p[0]], Extent::Revolved).is_err());
+ }
}
diff --git a/crates/model/src/sections.rs b/crates/model/src/sections.rs
index 3afc2f1..58a89d4 100644
--- a/crates/model/src/sections.rs
+++ b/crates/model/src/sections.rs
@@ -120,29 +120,90 @@ pub enum SectionKind {
Beam { profile: BeamProfile, orientation: [f64; 3] },
/// Axial members (rods, trusses) with a cross-sectional area.
Rod { area: f64 },
+ /// 2D continuum in plane stress (σzz = 0: thin plates loaded in their plane), of the given
+ /// thickness. The elements lie in the XY plane.
+ PlaneStress { thickness: f64 },
+ /// 2D continuum in plane strain (εzz = 0: a slice of a long body), of the given thickness
+ /// (1 by default, so results are per unit length). The elements lie in the XY plane.
+ PlaneStrain {
+ #[serde(default = "unit_thickness")]
+ thickness: f64,
+ },
+ /// Axisymmetric continuum: the elements lie in the XY half plane x ≥ 0, which is the r–z
+ /// plane of a body of revolution about the Y axis (x = radius, y = axial position; the hoop
+ /// direction is the third). Forces, reactions and masses are totals over the full
+ /// circumference (2π radians).
+ Axisymmetric {},
}
fn default_shell_points() -> u32 {
5
}
+fn unit_thickness() -> f64 {
+ 1.0
+}
+
+/// How a 2D continuum section extends out of its plane: what turns a length along an element
+/// edge into an area, and an element's area into a volume.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub enum Extent {
+ /// Plane stress and plane strain: a uniform thickness.
+ Thickness(f64),
+ /// Axisymmetric: the full circumference 2πr at radius r = x.
+ Revolved,
+}
+
+impl Extent {
+ /// The out-of-plane measure at a point: the thickness, or 2π x.
+ pub fn at(&self, xyz: [f64; 3]) -> f64 {
+ match *self {
+ Extent::Thickness(t) => t,
+ Extent::Revolved => 2.0 * std::f64::consts::PI * xyz[0],
+ }
+ }
+}
+
impl SectionKind {
/// The element dimension this section applies to.
pub fn dimension(&self) -> u8 {
match self {
SectionKind::Solid {} => 3,
- SectionKind::Shell { .. } => 2,
+ SectionKind::Shell { .. } | SectionKind::PlaneStress { .. } | SectionKind::PlaneStrain { .. } | SectionKind::Axisymmetric {} => 2,
SectionKind::Beam { .. } | SectionKind::Rod { .. } => 1,
}
}
- /// `solid`, `shell`, `beam` or `rod`.
+ /// Whether this is a continuum section (solid or 2D), the kinds with incompatible-mode
+ /// formulations.
+ pub fn is_continuum(&self) -> bool {
+ self.dimension() == 3 || self.is_planar()
+ }
+
+ /// Whether this is a 2D continuum section (plane stress, plane strain or axisymmetric).
+ pub fn is_planar(&self) -> bool {
+ self.extent().is_some()
+ }
+
+ /// How a 2D continuum section extends out of its plane (`None` for other kinds).
+ pub fn extent(&self) -> Option {
+ match *self {
+ SectionKind::PlaneStress { thickness } | SectionKind::PlaneStrain { thickness } => Some(Extent::Thickness(thickness)),
+ SectionKind::Axisymmetric {} => Some(Extent::Revolved),
+ _ => None,
+ }
+ }
+
+ /// `solid`, `shell`, `beam`, `rod`, `plane_stress`, `plane_strain` or `axisymmetric`.
pub fn name(&self) -> &'static str {
match self {
SectionKind::Solid {} => "solid",
SectionKind::Shell { .. } => "shell",
SectionKind::Beam { .. } => "beam",
SectionKind::Rod { .. } => "rod",
+ SectionKind::PlaneStress { .. } => "plane_stress",
+ SectionKind::PlaneStrain { .. } => "plane_strain",
+ SectionKind::Axisymmetric {} => "axisymmetric",
}
}
@@ -172,6 +233,10 @@ impl SectionKind {
SectionKind::Rod { area } => {
positive("area", *area)?;
}
+ SectionKind::PlaneStress { thickness } | SectionKind::PlaneStrain { thickness } => {
+ positive("thickness", *thickness)?;
+ }
+ SectionKind::Axisymmetric {} => {}
}
Ok(())
}
@@ -194,7 +259,8 @@ pub enum Integration {
pub struct Formulation {
#[serde(default)]
pub integration: Integration,
- /// Enhanced (incompatible) bending modes: linear hexahedra only.
+ /// Enhanced (incompatible) bending modes: linear hexahedra (solid sections) and linear
+ /// quadrilaterals (2D sections) only.
#[serde(default)]
pub incompatible_modes: bool,
/// Mixed displacement–pressure formulation for nearly incompressible materials: solids only.
@@ -209,8 +275,11 @@ impl Formulation {
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())));
+ if !kind.is_continuum() && self.incompatible_modes {
+ return Err(Error::Invalid(format!("incompatible modes are for solid and 2D sections, not {}", kind.name())));
+ }
+ if kind.dimension() != 3 && self.hybrid {
+ return Err(Error::Invalid(format!("hybrid formulations are for solid sections, not {}", kind.name())));
}
Ok(())
}
@@ -303,4 +372,26 @@ mod tests {
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());
}
+
+ #[test]
+ fn planar_sections() {
+ let strain: SectionKind = serde_json::from_value(serde_json::json!({"type": "plane_strain"})).unwrap();
+ assert_eq!(strain, SectionKind::PlaneStrain { thickness: 1.0 });
+ let stress: SectionKind = serde_json::from_value(serde_json::json!({"type": "plane_stress", "thickness": 2.5})).unwrap();
+ assert_eq!(stress.extent(), Some(Extent::Thickness(2.5)));
+ assert!(serde_json::from_value::(serde_json::json!({"type": "plane_stress"})).is_err(), "plane stress needs a thickness");
+ let axi: SectionKind = serde_json::from_value(serde_json::json!({"type": "axisymmetric"})).unwrap();
+ assert_eq!(axi.extent(), Some(Extent::Revolved));
+ assert!((Extent::Revolved.at([2.0, 7.0, 0.0]) - 4.0 * std::f64::consts::PI).abs() < 1e-15);
+ assert!(axi.is_planar() && stress.is_planar() && !SectionKind::Solid {}.is_planar());
+ assert!(!SectionKind::Shell { thickness: 1.0, offset: 0.0, points: 5 }.is_planar());
+ for bad in [SectionKind::PlaneStress { thickness: 0.0 }, SectionKind::PlaneStrain { thickness: f64::NAN }] {
+ assert!(bad.validate().is_err(), "{bad:?}");
+ }
+ let im = Formulation { incompatible_modes: true, ..Formulation::default() };
+ im.validate_for(&stress).unwrap();
+ im.validate_for(&axi).unwrap();
+ assert!(Formulation { hybrid: true, ..Formulation::default() }.validate_for(&strain).is_err());
+ assert!(im.validate_for(&SectionKind::Shell { thickness: 1.0, offset: 0.0, points: 5 }).is_err());
+ }
}
diff --git a/crates/model/src/sets.rs b/crates/model/src/sets.rs
index fff5207..4ee451f 100644
--- a/crates/model/src/sets.rs
+++ b/crates/model/src/sets.rs
@@ -10,10 +10,11 @@
//! 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;
+//! - an edge: the element edges on it (what 2D models load), their nodes, 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.
+//! those faces and their nodes; element edges: those edges and their nodes.
use std::collections::BTreeMap;
@@ -22,7 +23,8 @@ 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,
+ ElementEdge, 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.
@@ -40,6 +42,8 @@ pub struct Set {
pub elements: Vec,
#[serde(default, skip_serializing_if = "Vec::is_empty")]
pub faces: Vec,
+ #[serde(default, skip_serializing_if = "Vec::is_empty")]
+ pub edges: Vec,
}
impl Set {
@@ -49,35 +53,40 @@ impl Set {
}
pub fn is_empty(&self) -> bool {
- self.geometry.is_empty() && self.nodes.is_empty() && self.elements.is_empty() && self.faces.is_empty()
+ self.geometry.is_empty() && self.nodes.is_empty() && self.elements.is_empty() && self.faces.is_empty() && self.edges.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()));
+ return Err(Error::Invalid("a set needs geometry, nodes, elements, element faces or element edges".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 {
+ if self.nodes.len() > MAX_NODES || self.elements.len() > MAX_ELEMENTS || self.faces.len() > MAX_ELEMENTS || self.edges.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)) {
+ for &l in self.elements.iter().chain(self.faces.iter().map(|f| &f.element)).chain(self.edges.iter().map(|e| &e.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)));
}
+ if let Some(e) = self.edges.iter().find(|e| e.edge == 0) {
+ return Err(Error::Invalid(format!("element {} edge 0: edges are numbered from 1", e.element)));
+ }
for v in [&mut self.nodes, &mut self.elements] {
v.sort_unstable();
v.dedup();
}
self.faces.sort_unstable();
self.faces.dedup();
+ self.edges.sort_unstable();
+ self.edges.dedup();
self.geometry.sort();
self.geometry.dedup();
Ok(self)
@@ -134,6 +143,7 @@ impl Model {
let mut nodes = Vec::new();
let mut elements = Vec::new();
let mut faces = Vec::new();
+ let mut edges = Vec::new();
let mut missing = Vec::new();
// Geometry: mark the selected tags by role, then make one pass over the tagged elements.
@@ -192,8 +202,7 @@ impl Model {
}
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)));
+ edges.push(ElementEdge { element: label, edge });
}
}
}
@@ -224,16 +233,28 @@ impl Model {
None => missing.push(format!("element {} is not in the mesh", f.element)),
}
}
+ for &g in &set.edges {
+ match self.elements.get(&g.element) {
+ Some(e) if usize::from(g.edge) <= e.kind.edges().len() => edges.push(g),
+ Some(e) => missing.push(format!("element {} ({}) has no edge {}", g.element, e.kind.name(), g.edge)),
+ None => missing.push(format!("element {} is not in the mesh", g.element)),
+ }
+ }
elements.sort_unstable();
elements.dedup();
faces.sort_unstable();
faces.dedup();
+ edges.sort_unstable();
+ edges.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)));
}
}
+ for g in &edges {
+ nodes.extend(self.edge_nodes(*g).unwrap_or_default());
+ }
nodes.sort_unstable();
nodes.dedup();
// Long lists of missing labels are summarised.
@@ -242,7 +263,15 @@ impl Model {
missing.truncate(20);
missing.push(format!("and {more} more"));
}
- (ResolvedSet { nodes, elements, faces }, missing)
+ (ResolvedSet { nodes, elements, faces, edges }, missing)
+ }
+
+ /// The node labels of an element edge, in the edge's order ([`crate::ElementKind::edges`]:
+ /// its ends, then the mid-side node of a quadratic element).
+ pub fn edge_nodes(&self, edge: ElementEdge) -> Option> {
+ let e = self.elements.get(&edge.element)?;
+ let local = e.kind.edges().get(usize::from(edge.edge).checked_sub(1)?)?;
+ local.iter().map(|&i| e.nodes.get(i).copied()).collect()
}
/// The node labels of an element face, in the face's order ([`crate::ElementKind::faces`]).
diff --git a/crates/model/tests/definition.rs b/crates/model/tests/definition.rs
index 212bb6c..596815d 100644
--- a/crates/model/tests/definition.rs
+++ b/crates/model/tests/definition.rs
@@ -396,3 +396,113 @@ fn reference_point_labels_are_stable() {
"41 is a node of the remeshed bar"
);
}
+
+/// Two quad4 in the XY plane (x from `x0`, 2 × 1), tagged as face `area` of body `Plate` with the
+/// edges `left` (x = x0) and `right`, and a static step.
+fn plate(kind: SectionKind, x0: f64) -> Model {
+ let mut m = Model::default();
+ for (i, x) in [x0, x0 + 1.0, x0 + 2.0].iter().enumerate() {
+ m.set_node(i as u32 + 1, [*x, 0.0, 0.0]).unwrap();
+ m.set_node(i as u32 + 4, [*x, 1.0, 0.0]).unwrap();
+ }
+ m.set_element(1, ElementKind::Quad4, vec![1, 2, 5, 4]).unwrap();
+ m.set_element(2, ElementKind::Quad4, vec![2, 3, 6, 5]).unwrap();
+ let area = m.geometry_tag(&GeometryRef::face("Plate", "area").unwrap()).unwrap();
+ let left = m.geometry_tag(&GeometryRef::new("Plate", Entity::Edge("left".into())).unwrap()).unwrap();
+ let right = m.geometry_tag(&GeometryRef::new("Plate", Entity::Edge("right".into())).unwrap()).unwrap();
+ m.tag_element(1, area).unwrap();
+ m.tag_element(2, area).unwrap();
+ m.tag_element_edge(1, 4, left).unwrap();
+ m.tag_element_edge(2, 2, right).unwrap();
+ m.add_set("all", Set::of_geometry(vec![GeometryRef::face("Plate", "area").unwrap()])).unwrap();
+ m.add_set("left", Set::of_geometry(vec![GeometryRef::new("Plate", Entity::Edge("left".into())).unwrap()])).unwrap();
+ m.add_set("right", Set::of_geometry(vec![GeometryRef::new("Plate", Entity::Edge("right".into())).unwrap()])).unwrap();
+ m.add_material("Steel", Material { elastic: Some(Elastic::Isotropic { young: 2.1e5, poisson: 0.3 }), density: Some(7.85e-9), fluid: None })
+ .unwrap();
+ m.add_section("plate", Section { set: "all".into(), material: "Steel".into(), kind, formulation: Formulation::default() }).unwrap();
+ m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("left".into()) }).unwrap();
+ m.add_load("pull", Load::Traction { set: "right".into(), vector: [10.0, 0.0, 0.0] }).unwrap();
+ m.add_load("total", Load::Force { target: Target::Set("right".into()), vector: [0.0, 6.0, 0.0], distribution: Distribution::Total }).unwrap();
+ let mut s = Step::new("Pull", Procedure::Static {});
+ s.boundary_conditions.insert("clamp".into());
+ s.loads.extend(["pull".to_string(), "total".to_string()]);
+ m.add_step(s, None).unwrap();
+ m
+}
+
+#[test]
+fn two_dimensional_models_resolve_edges_and_check_their_rules() {
+ // Edge sets resolve to element edges; tractions spread over length × thickness.
+ let m = plate(SectionKind::PlaneStress { thickness: 2.0 }, 0.0);
+ let r = m.resolve().unwrap();
+ assert!(r.is_planar());
+ assert_eq!(r.planar_edges("right"), vec![ElementEdge { element: 2, edge: 2 }]);
+ assert_eq!(r.sets["left"].nodes, vec![1, 4]);
+ let pull: f64 = r.nodal_loads["pull"].iter().map(|(_, v)| v[0]).sum();
+ assert!((pull - 10.0 * 1.0 * 2.0).abs() < 1e-12, "{pull}");
+ assert_eq!(r.nodal_loads["total"], vec![(3, [0.0, 3.0, 0.0]), (6, [0.0, 3.0, 0.0])]);
+ // Axisymmetric: the right edge at r = 2 is a cylinder of area 2π·2·1.
+ let r = plate(SectionKind::Axisymmetric {}, 0.0).resolve().unwrap();
+ let pull: f64 = r.nodal_loads["pull"].iter().map(|(_, v)| v[0]).sum();
+ assert!((pull - 10.0 * 4.0 * std::f64::consts::PI).abs() < 1e-9, "{pull}");
+
+ let errors = |m: &Model| -> Vec { m.check().into_iter().filter(|d| d.severity == Severity::Error).map(|d| d.to_string()).collect() };
+ assert!(errors(&plate(SectionKind::Axisymmetric {}, -1.0)).iter().any(|e| e.contains("negative x")));
+ let mut m = plate(SectionKind::PlaneStrain { thickness: 1.0 }, 0.0);
+ m.set_node(6, [2.0, 1.0, 0.5]).unwrap();
+ m.set_element(1, ElementKind::Quad4, vec![1, 4, 5, 2]).unwrap();
+ let area = m.geometry_tag(&GeometryRef::face("Plate", "area").unwrap()).unwrap();
+ m.tag_element(1, area).unwrap();
+ m.add_load("lift", Load::Force { target: Target::Set("right".into()), vector: [0.0, 0.0, 1.0], distribution: Distribution::Total }).unwrap();
+ m.add_load("spin", Load::Centrifugal { angular_velocity: 10.0, point: [0.0; 3], axis: [0.0, 1.0, 0.0], set: None }).unwrap();
+ let e = errors(&m);
+ assert!(e.iter().any(|e| e.contains("off the XY plane")), "{e:#?}");
+ assert!(e.iter().any(|e| e.contains("clockwise")), "{e:#?}");
+ assert!(e.iter().any(|e| e.contains("load `lift`") && e.contains("Z component")), "{e:#?}");
+ assert!(e.iter().any(|e| e.contains("load `spin`") && e.contains("about an axis along Z")), "{e:#?}");
+ // A 2D section with a solid one; quad4 takes incompatible modes in a 2D section.
+ let mut m = plate(SectionKind::PlaneStress { thickness: 1.0 }, 0.0);
+ m.replace_section(
+ "plate",
+ Section {
+ set: "all".into(),
+ material: "Steel".into(),
+ kind: SectionKind::PlaneStress { thickness: 1.0 },
+ formulation: Formulation { incompatible_modes: true, ..Formulation::default() },
+ },
+ )
+ .unwrap();
+ assert!(errors(&m).is_empty(), "{:#?}", errors(&m));
+ for v in 1..=4 {
+ m.set_node(10 + v, [f64::from(v & 1), f64::from(v >> 1 & 1), 5.0]).unwrap();
+ }
+ m.set_element(9, ElementKind::Tet4, vec![11, 12, 13, 14]).unwrap();
+ m.add_set("tet", Set { elements: vec![9], ..Set::default() }).unwrap();
+ m.add_section("tet", Section { set: "tet".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() })
+ .unwrap();
+ assert!(errors(&m).iter().any(|e| e.contains("cannot be mixed")), "{:#?}", errors(&m));
+ // Pressure on an edge set of a 2D model is fine; on nothing it is an error.
+ let mut m = plate(SectionKind::PlaneStress { thickness: 1.0 }, 0.0);
+ m.add_load("press", Load::Pressure { set: "right".into(), magnitude: 1.0 }).unwrap();
+ m.add_load("bad", Load::Pressure { set: "all".into(), magnitude: 1.0 }).unwrap();
+ let e = errors(&m);
+ assert!(!e.iter().any(|e| e.contains("load `press`")) && e.iter().any(|e| e.contains("load `bad`")), "{e:#?}");
+ // Explicit element edges round-trip through JSON; edge 0 is an error.
+ let set = Set { edges: vec![ElementEdge { element: 2, edge: 2 }], ..Set::default() };
+ let json = serde_json::to_value(&set).unwrap();
+ assert_eq!(json, serde_json::json!({"edges": [{"element": 2, "edge": 2}]}));
+ assert!(m.add_set("zero", Set { edges: vec![ElementEdge { element: 2, edge: 0 }], ..Set::default() }).is_err());
+ assert!(m.add_set("five", Set { edges: vec![ElementEdge { element: 2, edge: 5 }], ..Set::default() }).is_ok());
+ assert!(m.resolve_set("five").is_err());
+ let json = serde_json::to_string(&m).unwrap();
+ assert_eq!(serde_json::from_str::(&json).unwrap(), m);
+ // Hostile centrifugal loads.
+ for bad in [
+ Load::Centrifugal { angular_velocity: 0.0, point: [0.0; 3], axis: [0.0, 1.0, 0.0], set: None },
+ Load::Centrifugal { angular_velocity: 1e200, point: [0.0; 3], axis: [0.0, 1.0, 0.0], set: None },
+ Load::Centrifugal { angular_velocity: 1.0, point: [f64::NAN, 0.0, 0.0], axis: [0.0, 1.0, 0.0], set: None },
+ Load::Centrifugal { angular_velocity: 1.0, point: [0.0; 3], axis: [0.0; 3], set: None },
+ ] {
+ assert!(m.add_load("hostile", bad.clone()).is_err(), "{bad:?}");
+ }
+}
diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs
index a143e59..700c0f4 100644
--- a/crates/solver/src/analysis.rs
+++ b/crates/solver/src/analysis.rs
@@ -1,5 +1,6 @@
//! Linear static and frequency analysis of a model's step, and their results.
+use modelcraft_fem::ElementType;
use modelcraft_fem::element::{Extrapolation, extrapolation};
use modelcraft_fem::quadrature::Rule;
use modelcraft_fem::{DofKind, DofLayout, PointStress, von_mises};
@@ -134,14 +135,15 @@ fn assemble(mesh: &Mesh, graph: &NodeGraph, mass: Option) -> Result add_scalar(mm, layout, &e.nodes, solid.dofs().node, &solid.mass(rho).map_err(at)?.data)?,
- MassKind::Lumped => add_diagonal(mm, layout, &e.nodes, solid.dofs().node, &solid.lumped_mass(rho).map_err(at)?)?,
+ MassKind::Consistent => add_scalar(mm, layout, &e.nodes, kinds, &kernel.mass(rho).map_err(at)?.data)?,
+ MassKind::Lumped => add_diagonal(mm, layout, &e.nodes, kinds, &kernel.lumped_mass(rho).map_err(at)?)?,
}
}
}
@@ -243,22 +245,22 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul
.map(|(p, (d, r))| PointResult { name: p.name.clone(), label: p.label, displacement: *d, reaction: *r })
.collect();
// Stresses: one extrapolation per element type and rule.
- let mut extrapolations: Vec<(ElementKind, Rule, Extrapolation)> = Vec::new();
+ let mut extrapolations: Vec<(&'static ElementType, Rule, Extrapolation)> = Vec::new();
let mut elements = Vec::with_capacity(mesh.elements.len());
let mut sum = vec![[0.0; 6]; mesh.labels.len()];
let mut count = vec![0u32; mesh.labels.len()];
for e in &mesh.elements {
- let solid = e.solid()?;
+ let kernel = e.kernel()?;
let at = |err| Error::element(e.label, err);
let ue: Vec<[f64; 3]> = e.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect();
- let points = solid.stresses(&e.d, &ue).map_err(at)?;
+ let points = kernel.stresses(&e.d, &ue).map_err(at)?;
let values: Vec<[f64; 6]> = points.iter().map(|p| p.stress).collect();
- let rule = solid.rule();
- let cached = extrapolations.iter().position(|(k, r, _)| *k == e.kind && *r == rule);
+ let (ty, rule) = (kernel.ty(), kernel.rule());
+ let cached = extrapolations.iter().position(|(t, r, _)| std::ptr::eq(*t, ty) && *r == rule);
let index = match cached {
Some(i) => i,
None => {
- extrapolations.push((e.kind, rule, extrapolation(solid.ty, rule).map_err(at)?));
+ extrapolations.push((ty, rule, extrapolation(ty, rule).map_err(at)?));
extrapolations.len() - 1
}
};
@@ -341,7 +343,7 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt
let Assembled { k, m } = assemble(&mesh, &graph, Some(options.mass))?;
let m = m.ok_or_else(|| Error::Internal("mass matrix".into()))?;
// Every element's mass entries sum to its mass, once per direction.
- let total_mass = m.val.iter().sum::() / 3.0;
+ let total_mass = m.val.iter().sum::() / if mesh.theory.is_some() { 2.0 } else { 3.0 };
let neq = red.neq();
if neq == 0 {
return Err(Error::Invalid("every degree of freedom is held; there is nothing to vibrate".into()));
diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs
index 3d8d828..8dbac28 100644
--- a/crates/solver/src/system.rs
+++ b/crates/solver/src/system.rs
@@ -3,18 +3,20 @@
//! load vector of a step.
//!
//! Degrees of freedom ([`DofLayout`]): the kinds each mesh node's elements declare (three
-//! translations for solids), nodes by label, then six per coupled reference point (points by
+//! translations for solids, Ux and Uy for 2D elements), nodes by label, then six per coupled reference point (points by
//! name): translations along X, Y, Z and rotations about them.
use modelcraft_fem::element::{Formulation, Solid};
-use modelcraft_fem::{DofKind, DofLayout, DofSet, Elasticity, elasticity, face};
+use modelcraft_fem::plane::{Plane, Theory};
+use modelcraft_fem::{DofKind, DofLayout, DofSet, Elasticity, ElementDofs, PointStress, elasticity, face, plane};
+use modelcraft_model::nodal::edge_pressure;
use std::collections::BTreeSet;
-use modelcraft_model::{ElementKind, Load, Model, Resolved, SectionKind, Step, Target};
+use modelcraft_model::{ElementKind, Load, Model, Resolved, SectionKind, Step, Target, centrifugal_acceleration};
use crate::{Error, Result};
-/// A solid element ready for the kernels.
+/// An element ready for the kernels: a solid, or a 2D element (`theory`).
#[derive(Debug, Clone)]
pub struct Elem {
pub label: u32,
@@ -25,11 +27,93 @@ pub struct Elem {
pub formulation: Formulation,
pub d: Elasticity,
pub density: Option,
+ /// The 2D theory of an element of a 2D section; `None` for solids.
+ pub theory: Option,
+}
+
+/// An element's kernel: a solid or a 2D element.
+#[derive(Debug, Clone, Copy)]
+pub enum Kernel<'a> {
+ Solid(Solid<'a>),
+ Plane(Plane<'a>),
+}
+
+impl Kernel<'_> {
+ /// The element's unknowns, in the order of its matrices.
+ pub fn dofs(&self) -> ElementDofs {
+ match self {
+ Kernel::Solid(s) => s.dofs(),
+ Kernel::Plane(p) => p.dofs(),
+ }
+ }
+
+ pub fn stiffness(&self, d: &Elasticity) -> modelcraft_fem::Result {
+ match self {
+ Kernel::Solid(s) => s.stiffness(d),
+ Kernel::Plane(p) => p.stiffness(d),
+ }
+ }
+
+ pub fn mass(&self, density: f64) -> modelcraft_fem::Result {
+ match self {
+ Kernel::Solid(s) => s.mass(density),
+ Kernel::Plane(p) => p.mass(density),
+ }
+ }
+
+ pub fn lumped_mass(&self, density: f64) -> modelcraft_fem::Result> {
+ match self {
+ Kernel::Solid(s) => s.lumped_mass(density),
+ Kernel::Plane(p) => p.lumped_mass(density),
+ }
+ }
+
+ /// Consistent nodal forces of a uniform force per volume, as 3-vectors (Z zero for 2D).
+ pub fn body_force(&self, b: [f64; 3]) -> modelcraft_fem::Result> {
+ match self {
+ Kernel::Solid(s) => s.body_force(b),
+ Kernel::Plane(p) => Ok(p.body_force_field(|_| b)?.into_iter().map(|v| [v[0], v[1], 0.0]).collect()),
+ }
+ }
+
+ /// Consistent nodal forces of a force per volume b(x), as 3-vectors (Z zero for 2D).
+ pub fn body_force_field(&self, b: impl Fn([f64; 3]) -> [f64; 3]) -> modelcraft_fem::Result> {
+ match self {
+ Kernel::Solid(s) => s.body_force_field(b),
+ Kernel::Plane(p) => Ok(p.body_force_field(b)?.into_iter().map(|v| [v[0], v[1], 0.0]).collect()),
+ }
+ }
+
+ /// Integration-point stresses for nodal displacements given as 3-vectors.
+ pub fn stresses(&self, d: &Elasticity, u: &[[f64; 3]]) -> modelcraft_fem::Result> {
+ match self {
+ Kernel::Solid(s) => s.stresses(d, u),
+ Kernel::Plane(p) => p.stresses(d, &u.iter().map(|v| [v[0], v[1]]).collect::>()),
+ }
+ }
+
+ pub fn rule(&self) -> modelcraft_fem::quadrature::Rule {
+ match self {
+ Kernel::Solid(s) => s.rule(),
+ Kernel::Plane(p) => p.rule(),
+ }
+ }
+
+ pub fn ty(&self) -> &'static modelcraft_fem::ElementType {
+ match self {
+ Kernel::Solid(s) => s.ty,
+ Kernel::Plane(p) => p.ty,
+ }
+ }
}
impl Elem {
- pub fn solid(&self) -> Result> {
- Solid::new(self.kind, &self.xyz, self.formulation).map_err(|e| Error::element(self.label, e))
+ pub fn kernel(&self) -> Result> {
+ let k = match self.theory {
+ None => Solid::new(self.kind, &self.xyz, self.formulation).map(Kernel::Solid),
+ Some(t) => Plane::new(self.kind, &self.xyz, self.formulation, t).map(Kernel::Plane),
+ };
+ k.map_err(|e| Error::element(self.label, e))
}
}
@@ -52,6 +136,8 @@ pub struct Mesh {
pub elements: Vec,
pub points: Vec,
pub layout: DofLayout,
+ /// The 2D theory when this is a 2D model (a model is 2D as a whole).
+ pub theory: Option,
}
impl Mesh {
@@ -102,13 +188,19 @@ impl Mesh {
/// Build it from the model and its resolved view, refusing what the solver cannot represent.
pub fn new(model: &Model, resolved: &Resolved, max_dofs: usize) -> Result {
let mut elements = Vec::new();
+ // A model is 2D as a whole (the model check refuses a mix).
+ let planar = model.sections().values().find_map(|s| Theory::from_section(&s.kind));
for (name, section) in model.sections() {
- if !matches!(section.kind, SectionKind::Solid {}) {
+ let theory = Theory::from_section(§ion.kind);
+ if !matches!(section.kind, SectionKind::Solid {}) && theory.is_none() {
return Err(Error::NotSupported(format!(
- "{} sections (section `{name}`); the solver has solid elements so far",
+ "{} sections (section `{name}`); the solver has solid and 2D (plane stress, plane strain, axisymmetric) elements so far",
section.kind.name()
)));
}
+ if theory.is_some() != planar.is_some() {
+ return Err(Error::NotSupported("2D and solid sections in one model".into()));
+ }
let material = model
.material(§ion.material)
.ok_or_else(|| Error::Invalid(format!("section `{name}`: material `{}` is missing", section.material)))?;
@@ -118,10 +210,17 @@ impl Mesh {
let set = resolved.sets.get(§ion.set).ok_or_else(|| Error::Invalid(format!("set `{}` is not resolved", section.set)))?;
for &label in &set.elements {
let el = model.element(label).ok_or_else(|| Error::Invalid(format!("element {label} is missing")))?;
- let formulation = Formulation::from_section(el.kind, §ion.formulation).map_err(|e| Error::element(label, e))?;
- elements.push((label, el.kind, el.nodes.clone(), formulation, d, material.density));
+ let (formulation, dofs) = match theory {
+ None => Formulation::from_section(el.kind, §ion.formulation).map(|f| (f, f.dofs())),
+ Some(_) => plane::formulation(el.kind, §ion.formulation).map(|f| (f, ElementDofs::new(DofSet::PLANE))),
+ }
+ .map_err(|e| Error::element(label, e))?;
+ elements.push((label, el.kind, el.nodes.clone(), formulation, d, material.density, theory, dofs));
}
}
+ if let (Some(_), Some(name)) = (planar, model.couplings().keys().next()) {
+ return Err(Error::NotSupported(format!("couplings in 2D models (coupling `{name}`)")));
+ }
elements.sort_by_key(|e| e.0);
let mut labels: Vec = elements.iter().flat_map(|e| e.2.iter().copied()).collect();
labels.sort_unstable();
@@ -136,7 +235,7 @@ impl Mesh {
.collect::>()?;
let mut layout = DofLayout::builder();
for e in &elements {
- layout.element(&e.2, e.3.dofs());
+ layout.element(&e.2, e.7);
}
for p in &points {
layout.point(p.label, DofSet::STRUCTURAL);
@@ -150,11 +249,11 @@ impl Mesh {
.iter()
.map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| Error::Invalid(format!("node {l} is missing"))))
.collect::>()?;
- let mut mesh = Mesh { labels, xyz, elements: Vec::with_capacity(elements.len()), points, layout };
- for (label, kind, nodes, formulation, d, density) in elements {
+ let mut mesh = Mesh { labels, xyz, elements: Vec::with_capacity(elements.len()), points, layout, theory: planar };
+ for (label, kind, nodes, formulation, d, density, theory, _) in elements {
let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect();
let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect();
- mesh.elements.push(Elem { label, kind, nodes: idx, xyz: exyz, formulation, d, density });
+ mesh.elements.push(Elem { label, kind, nodes: idx, xyz: exyz, formulation, d, density, theory });
}
Ok(mesh)
}
@@ -205,11 +304,31 @@ pub fn prescribed(
}
}
}
+ // Axisymmetric: a node on the axis cannot move radially.
+ if mesh.theory == Some(Theory::Axisymmetric) {
+ for (n, p) in mesh.xyz.iter().enumerate() {
+ let index = mesh.dof(n, DofKind::Ux)?;
+ if p[0] != 0.0 || !active.get(index).copied().unwrap_or(false) {
+ continue;
+ }
+ let slot = out.get_mut(index).ok_or_else(|| Error::Internal("dof index".into()))?;
+ match *slot {
+ Some(v) if v != 0.0 && values => {
+ return Err(Error::Invalid(format!(
+ "{} is on the axis, where the radial displacement is zero, not {v}",
+ mesh.describe_dof(index)
+ )));
+ }
+ _ => *slot = Some(0.0),
+ }
+ }
+ }
Ok(out)
}
-/// Add a vector to the three degrees of freedom `dofs`.
-fn add(f: &mut [f64], dofs: [usize; 3], v: [f64; 3]) {
+/// Add a vector's components to the degrees of freedom `dofs` (the first `dofs.len()`
+/// components).
+fn add(f: &mut [f64], dofs: &[usize], v: [f64; 3]) {
for (d, c) in dofs.iter().zip(v) {
if let Some(x) = f.get_mut(*d) {
*x += c;
@@ -217,15 +336,20 @@ fn add(f: &mut [f64], dofs: [usize; 3], v: [f64; 3]) {
}
}
-/// The translations (or rotations) of layout slot `slot`, if it has them.
-fn vector_dofs(mesh: &Mesh, slot: usize, rotation: bool) -> Option<[usize; 3]> {
+/// The translations (or rotations) of layout slot `slot`, if it has them: X, Y, Z, or X and Y
+/// for a node of a 2D model.
+fn vector_dofs(mesh: &Mesh, slot: usize, rotation: bool) -> Option> {
let kind = |c| if rotation { DofKind::rotation(c) } else { DofKind::translation(c) };
- let [a, b, c] = [0, 1, 2].map(|c| kind(c).and_then(|k| mesh.layout.dof(slot, k)));
- Some([a?, b?, c?])
+ let dofs: Vec> = [0, 1, 2].map(|c| kind(c).and_then(|k| mesh.layout.dof(slot, k))).into();
+ match dofs.as_slice() {
+ [Some(a), Some(b), Some(c)] => Some(vec![*a, *b, *c]),
+ [Some(a), Some(b), None] if !rotation && mesh.theory.is_some() => Some(vec![*a, *b]),
+ _ => None,
+ }
}
/// The translations of mesh node `node`.
-fn translations(mesh: &Mesh, node: usize) -> Result<[usize; 3]> {
+fn translations(mesh: &Mesh, node: usize) -> Result> {
vector_dofs(mesh, node, false).ok_or_else(|| Error::Internal(format!("mesh node {node} has no translations")))
}
@@ -259,7 +383,28 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ
)));
}
}
- add(&mut f, dofs, *v);
+ add(&mut f, &dofs, *v);
+ }
+ }
+ Load::Pressure { set, magnitude } if mesh.theory.is_some() => {
+ if resolved.sets.get(set).is_some_and(|s| !s.faces.is_empty()) {
+ return Err(Error::Invalid(format!(
+ "load `{name}`: a 2D model's pressures act on element edges, and set `{set}` has element faces"
+ )));
+ }
+ for g in resolved.planar_edges(set) {
+ let e = mesh.element(g.element).ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} is missing", g.element)))?;
+ let labels = model
+ .edge_nodes(g)
+ .ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} edge {} is missing", g.element, g.edge)))?;
+ let idx: Vec = labels.iter().filter_map(|l| mesh.node(*l)).collect();
+ let xyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect();
+ let theory = e.theory.ok_or_else(|| Error::Internal("2D element without a theory".into()))?;
+ let forces = edge_pressure(&xyz, theory.extent(), *magnitude)
+ .map_err(|err| Error::element(g.element, modelcraft_fem::Error::Invalid(err.to_string())))?;
+ for (i, v) in idx.iter().zip(forces) {
+ add(&mut f, &translations(mesh, *i)?, v);
+ }
}
}
Load::Pressure { set, magnitude } => {
@@ -275,7 +420,7 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ
let xyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect();
let forces = face::pressure(&xyz, *magnitude).map_err(|e| Error::element(fc.element, e))?;
for (i, v) in idx.iter().zip(forces) {
- add(&mut f, translations(mesh, *i)?, v);
+ add(&mut f, &translations(mesh, *i)?, v);
}
}
}
@@ -295,15 +440,31 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ
body(&mut f, mesh, e, *vector)?;
}
}
+ Load::Centrifugal { angular_velocity, point, axis, set } => {
+ let chosen: Option<&[u32]> = set.as_ref().map(|s| resolved.sets.get(s).map(|r| r.elements.as_slice()).unwrap_or_default());
+ for e in &mesh.elements {
+ if chosen.is_some_and(|c| c.binary_search(&e.label).is_err()) {
+ continue;
+ }
+ let rho = e.density.ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} has no density", e.label)))?;
+ let forces = e
+ .kernel()?
+ .body_force_field(|x| centrifugal_acceleration(*angular_velocity, *point, *axis, x).map(|a| rho * a))
+ .map_err(|err| Error::element(e.label, err))?;
+ for (node, v) in e.nodes.iter().zip(forces) {
+ add(&mut f, &translations(mesh, *node)?, v);
+ }
+ }
+ }
}
}
Ok(f)
}
fn body(f: &mut [f64], mesh: &Mesh, e: &Elem, b: [f64; 3]) -> Result<()> {
- let forces = e.solid()?.body_force(b).map_err(|err| Error::element(e.label, err))?;
+ let forces = e.kernel()?.body_force(b).map_err(|err| Error::element(e.label, err))?;
for (node, v) in e.nodes.iter().zip(forces) {
- add(f, translations(mesh, *node)?, v);
+ add(f, &translations(mesh, *node)?, v);
}
Ok(())
}
diff --git a/crates/solver/tests/plane.rs b/crates/solver/tests/plane.rs
new file mode 100644
index 0000000..c6b0121
--- /dev/null
+++ b/crates/solver/tests/plane.rs
@@ -0,0 +1,567 @@
+//! The analytic verification suite of the 2D elements (plane stress, plane strain,
+//! axisymmetric): patch tests, a traction patch, a cantilever, Kirsch's plate, Lamé's cylinder
+//! and rotating discs and cylinders against closed-form solutions. Run with `--nocapture` to see
+//! the errors.
+
+use std::collections::BTreeMap;
+use std::f64::consts::PI;
+
+use modelcraft_model::{
+ BoundaryCondition, Elastic, ElementEdge, ElementKind, Formulation, Integration, Load, Material, Model, Procedure, Section, SectionKind, Set,
+ Step, Symmetry, Target,
+};
+use modelcraft_solver::{FrequencyOptions, Limits, StaticResult, solve_frequency, solve_static};
+
+/// A structured 2D mesh of `n` cells on [0, 1]², mapped by `map`; nodes on the doubled lattice so
+/// quadratic elements take their mid-side nodes from the map.
+struct Grid {
+ model: Model,
+ param: BTreeMap,
+}
+
+fn near(a: f64, b: f64) -> bool {
+ (a - b).abs() < 1e-9
+}
+
+fn grid(n: [usize; 2], kind: ElementKind, map: impl Fn([f64; 2]) -> [f64; 2]) -> Grid {
+ let size = [2 * n[0] + 1, 2 * n[1] + 1];
+ let mut model = Model::default();
+ let mut param = BTreeMap::new();
+ let mut node = |model: &mut Model, i: [usize; 2]| -> u32 {
+ let l = (1 + i[0] + size[0] * i[1]) as u32;
+ param.entry(l).or_insert_with(|| {
+ let p = [i[0] as f64 / (2 * n[0]) as f64, i[1] as f64 / (2 * n[1]) as f64];
+ let x = map(p);
+ model.set_node(l, [x[0], x[1], 0.0]).unwrap();
+ p
+ });
+ l
+ };
+ let mut next = 1u32;
+ for cy in 0..n[1] {
+ for cx in 0..n[0] {
+ let at = |d: [usize; 2]| [2 * cx + d[0], 2 * cy + d[1]];
+ let polygons: Vec> = match kind {
+ ElementKind::Quad4 | ElementKind::Quad8 => vec![vec![[0, 0], [2, 0], [2, 2], [0, 2]]],
+ _ => vec![vec![[0, 0], [2, 0], [2, 2]], vec![[0, 0], [2, 2], [0, 2]]],
+ };
+ for corners in polygons {
+ let mut nodes: Vec = corners.iter().map(|c| node(&mut model, at(*c))).collect();
+ if kind.is_quadratic() {
+ for e in kind.edges() {
+ let (a, b) = (corners[e[0]], corners[e[1]]);
+ nodes.push(node(&mut model, at([(a[0] + b[0]) / 2, (a[1] + b[1]) / 2])));
+ }
+ }
+ model.set_element(next, kind, nodes).unwrap();
+ next += 1;
+ }
+ }
+ }
+ Grid { model, param }
+}
+
+impl Grid {
+ fn nodes_where(&self, f: impl Fn([f64; 2]) -> bool) -> Vec {
+ self.param.iter().filter(|(_, p)| f(**p)).map(|(l, _)| *l).collect()
+ }
+
+ fn edges_where(&self, f: impl Fn([f64; 2]) -> bool) -> Vec {
+ let mut out = Vec::new();
+ for (label, e) in self.model.elements() {
+ for (i, edge) in e.kind.edges().iter().enumerate() {
+ if edge.iter().all(|k| f(self.param[&e.nodes[*k]])) {
+ out.push(ElementEdge { element: *label, edge: i as u8 + 1 });
+ }
+ }
+ }
+ out
+ }
+
+ fn define(&mut self, kind: SectionKind, young: f64, nu: f64, density: Option, formulation: Formulation, procedure: Procedure) {
+ let all: Vec = self.model.elements().keys().copied().collect();
+ self.model.add_set("all", Set { elements: all, ..Set::default() }).unwrap();
+ self.model.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young, poisson: nu }), density, fluid: None }).unwrap();
+ self.model.add_section("plane", Section { set: "all".into(), material: "Mat".into(), kind, formulation }).unwrap();
+ self.model.add_step(Step::new("Step", procedure), None).unwrap();
+ }
+
+ fn bc(&mut self, name: &str, bc: BoundaryCondition) {
+ self.model.add_boundary_condition(name, bc).unwrap();
+ self.model.set_boundary_condition_active("Step", name, true).unwrap();
+ }
+
+ fn load(&mut self, name: &str, load: Load) {
+ self.model.add_load(name, load).unwrap();
+ self.model.set_load_active("Step", name, true).unwrap();
+ }
+
+ fn symmetry(&mut self, name: &str, plane: Symmetry, f: impl Fn([f64; 2]) -> bool) {
+ let nodes = self.nodes_where(f);
+ self.model.add_set(name, Set { nodes, ..Set::default() }).unwrap();
+ self.bc(name, BoundaryCondition::Symmetry { target: Target::Set(name.into()), plane });
+ }
+
+ fn fixed(&mut self, name: &str, f: impl Fn([f64; 2]) -> bool) {
+ let nodes = self.nodes_where(f);
+ self.model.add_set(name, Set { nodes, ..Set::default() }).unwrap();
+ self.bc(name, BoundaryCondition::Fixed { target: Target::Set(name.into()) });
+ }
+
+ fn edge_set(&mut self, name: &str, f: impl Fn([f64; 2]) -> bool) {
+ let edges = self.edges_where(f);
+ self.model.add_set(name, Set { edges, ..Set::default() }).unwrap();
+ }
+
+ fn solve(&self) -> StaticResult {
+ solve_static(&self.model, Some("Step"), &Limits::default()).unwrap()
+ }
+}
+
+fn at(r: &StaticResult, label: u32) -> usize {
+ r.nodes.binary_search(&label).unwrap()
+}
+
+fn full() -> Formulation {
+ Formulation::default()
+}
+
+fn reduced() -> Formulation {
+ Formulation { integration: Integration::Reduced, ..Formulation::default() }
+}
+
+fn incompatible() -> Formulation {
+ Formulation { incompatible_modes: true, ..Formulation::default() }
+}
+
+fn variants() -> Vec<(ElementKind, Formulation, &'static str)> {
+ vec![
+ (ElementKind::Tri3, full(), "tri3"),
+ (ElementKind::Tri6, full(), "tri6"),
+ (ElementKind::Quad4, full(), "quad4"),
+ (ElementKind::Quad4, incompatible(), "quad4 incompatible"),
+ (ElementKind::Quad8, full(), "quad8"),
+ (ElementKind::Quad8, reduced(), "quad8 reduced"),
+ ]
+}
+
+const STRESS: SectionKind = SectionKind::PlaneStress { thickness: 0.5 };
+const STRAIN: SectionKind = SectionKind::PlaneStrain { thickness: 1.0 };
+const AXI: SectionKind = SectionKind::Axisymmetric {};
+
+/// Displacement patch test on a distorted 2 × 2 patch over [1, 2] × [0, 1]: a linear
+/// displacement prescribed on the boundary is reproduced inside, with exactly constant strain
+/// and stress (plane: a general linear field; axisymmetric: u_r = a r, u_z = b z, the
+/// constant-strain state a body of revolution can be in without body forces; a constant shear
+/// γ_rz would need an axial body force τ/r).
+#[test]
+fn patch_test_every_element_and_theory() {
+ let (young, nu) = (1000.0, 0.25);
+ let map = |p: [f64; 2]| {
+ let b = (PI * p[0]).sin() * (PI * p[1]).sin();
+ [1.0 + p[0] + 0.15 * b, p[1] - 0.1 * b]
+ };
+ let lambda = young * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
+ let mu = young / (2.0 * (1.0 + nu));
+ for theory in [STRESS, STRAIN, AXI] {
+ // u = field(x); strain (xx, yy, zz, xy).
+ let (field, strain): (fn([f64; 3]) -> [f64; 2], [f64; 4]) = if theory == AXI {
+ (|x| [1e-3 * x[0], -5e-4 * x[1]], [1e-3, -5e-4, 1e-3, 0.0])
+ } else {
+ (|x| [1e-3 * x[0] + 2e-4 * x[1], 5e-4 * x[0] - 2e-3 * x[1]], [1e-3, -2e-3, 0.0, 7e-4])
+ };
+ let mut e = strain;
+ if let SectionKind::PlaneStress { .. } = theory {
+ e[2] = -nu / (1.0 - nu) * (e[0] + e[1]);
+ }
+ let tr = e[0] + e[1] + e[2];
+ let stress = [lambda * tr + 2.0 * mu * e[0], lambda * tr + 2.0 * mu * e[1], lambda * tr + 2.0 * mu * e[2], mu * e[3]];
+ for (kind, formulation, name) in variants() {
+ let mut g = grid([2, 2], kind, map);
+ if theory == AXI {
+ // Straight sides: on a curved quadratic element the axisymmetric integrand
+ // (r ∂N/∂r + N) det J is rational, so no Gauss rule integrates it exactly and the
+ // patch test holds only to the rule's accuracy.
+ straighten(&mut g.model);
+ }
+ g.define(theory, young, nu, None, formulation, Procedure::Static {});
+ let boundary = g.nodes_where(|p| near(p[0], 0.0) || near(p[0], 1.0) || near(p[1], 0.0) || near(p[1], 1.0));
+ for (i, label) in boundary.into_iter().enumerate() {
+ let u = field(g.model.node(label).unwrap().xyz);
+ let set = format!("b{i}");
+ g.model.add_set(&set, Set { nodes: vec![label], ..Set::default() }).unwrap();
+ g.bc(
+ &set,
+ BoundaryCondition::Displacement { target: Target::Set(set.clone()), values: [Some(u[0]), Some(u[1]), None, None, None, None] },
+ );
+ }
+ let r = g.solve();
+ let (mut eu, mut es): (f64, f64) = (0.0, 0.0);
+ for (label, u) in r.nodes.iter().zip(&r.displacement) {
+ let want = field(g.model.node(*label).unwrap().xyz);
+ eu = eu.max((u[0] - want[0]).abs().max((u[1] - want[1]).abs()));
+ assert_eq!(u[2], 0.0);
+ }
+ for el in &r.elements {
+ for p in &el.points {
+ for (got, want) in p.stress.iter().zip(&stress) {
+ es = es.max((got - want).abs() / young);
+ }
+ }
+ }
+ eprintln!("patch {} {name}: displacement error {eu:.1e}, stress error {es:.1e}", theory.name());
+ assert!(eu < 1e-15 && es < 1e-14, "{} {name}: {eu:e} {es:e}", theory.name());
+ }
+ }
+}
+
+/// Put every mid-side node at the middle of its edge's ends.
+fn straighten(model: &mut Model) {
+ let elements: Vec<(ElementKind, Vec)> = model.elements().values().map(|e| (e.kind, e.nodes.clone())).collect();
+ for (kind, nodes) in elements {
+ for e in kind.edges() {
+ if let [a, b, m] = **e {
+ let (pa, pb) = (model.node(nodes[a]).unwrap().xyz, model.node(nodes[b]).unwrap().xyz);
+ model.set_node(nodes[m], [0, 1, 2].map(|k| (pa[k] + pb[k]) / 2.0)).unwrap();
+ }
+ }
+ }
+}
+
+/// A uniform traction on the free edge of a plate (plane stress), a slab (plane strain) and a
+/// tube's end (axisymmetric), held by symmetry: the exact stress and displacement everywhere.
+#[test]
+fn uniform_traction_is_exact() {
+ let (young, nu, s0) = (1000.0, 0.3, 2.0);
+ for theory in [STRESS, STRAIN, AXI] {
+ for (kind, formulation, name) in variants() {
+ // Axisymmetric: a tube r ∈ [1, 2], pulled along its axis (Y); plane: pulled along X.
+ let mut g = grid([3, 2], kind, |p| [1.0 + p[0], 2.0 * p[1]]);
+ g.define(theory, young, nu, None, formulation, Procedure::Static {});
+ if theory == AXI {
+ g.symmetry("bottom", Symmetry::Y, |p| near(p[1], 0.0));
+ g.edge_set("top", |p| near(p[1], 1.0));
+ g.load("pull", Load::Traction { set: "top".into(), vector: [0.0, s0, 0.0] });
+ } else {
+ g.symmetry("left", Symmetry::X, |p| near(p[0], 0.0));
+ g.symmetry("bottom", Symmetry::Y, |p| near(p[1], 0.0));
+ g.edge_set("right", |p| near(p[0], 1.0));
+ g.load("pull", Load::Pressure { set: "right".into(), magnitude: -s0 });
+ }
+ let r = g.solve();
+ // Strains: along the load s/E' and across −ν' s/E'.
+ let (along, across) = match theory {
+ SectionKind::PlaneStrain { .. } => (s0 * (1.0 - nu * nu) / young, -nu * (1.0 + nu) * s0 / young),
+ _ => (s0 / young, -nu * s0 / young),
+ };
+ let mut eu: f64 = 0.0;
+ for (label, u) in r.nodes.iter().zip(&r.displacement) {
+ let x = g.model.node(*label).unwrap().xyz;
+ let want = if theory == AXI { [across * x[0], along * x[1]] } else { [along * (x[0] - 1.0), across * x[1]] };
+ eu = eu.max((u[0] - want[0]).abs().max((u[1] - want[1]).abs()));
+ }
+ let total: f64 = r.reaction_total.iter().map(|c| c.abs()).sum();
+ let area = if theory == AXI {
+ PI * (4.0 - 1.0)
+ } else {
+ 2.0 * if let SectionKind::PlaneStress { thickness } = theory { thickness } else { 1.0 }
+ };
+ eprintln!("traction {} {name}: displacement error {eu:.1e}, reaction {total:.6} vs {:.6}", theory.name(), s0 * area);
+ assert!(eu < 1e-13, "{} {name}: {eu:e}", theory.name());
+ assert!((total - s0 * area).abs() < 1e-11 * s0 * area, "{} {name}: {total} vs {}", theory.name(), s0 * area);
+ }
+ }
+}
+
+/// Cantilever L × h of thickness t in plane stress, clamped at x = 0, a shear traction on the free
+/// end; Timoshenko's deflection δ = P L³/(3 E I) + P L/(κ G A), κ = 5/6, with ν = 0.
+fn cantilever_tip(kind: ElementKind, formulation: Formulation, n: [usize; 2]) -> f64 {
+ let (l, h, t, young, p) = (10.0, 1.0, 1.0, 1000.0, 1.0);
+ let mut g = grid(n, kind, |q| [l * q[0], h * q[1]]);
+ g.define(SectionKind::PlaneStress { thickness: t }, young, 0.0, None, formulation, Procedure::Static {});
+ g.fixed("root", |q| near(q[0], 0.0));
+ g.edge_set("end", |q| near(q[0], 1.0));
+ g.load("shear", Load::Traction { set: "end".into(), vector: [0.0, -p / (h * t), 0.0] });
+ let r = g.solve();
+ let end = g.nodes_where(|q| near(q[0], 1.0));
+ -end.iter().map(|l| r.displacement[at(&r, *l)][1]).sum::() / end.len() as f64
+}
+
+#[test]
+fn cantilever_in_plane_stress() {
+ let (l, h, young, p) = (10.0f64, 1.0f64, 1000.0, 1.0);
+ let exact = p * l.powi(3) / (3.0 * young * h.powi(3) / 12.0) + p * l / (5.0 / 6.0 * young / 2.0 * h);
+ let err = |d: f64| (d - exact).abs() / exact;
+ let study = |kind, f, meshes: &[[usize; 2]], name: &str| -> Vec {
+ meshes
+ .iter()
+ .map(|n| {
+ let d = cantilever_tip(kind, f, *n);
+ eprintln!("cantilever {name} {n:?}: {d:.5} vs {exact:.5} (error {:.2}%)", 100.0 * err(d));
+ err(d)
+ })
+ .collect()
+ };
+ let q8 = study(ElementKind::Quad8, full(), &[[10, 1], [20, 2]], "quad8");
+ assert!(q8[1] < 0.005, "{q8:?}");
+ let q8r = study(ElementKind::Quad8, reduced(), &[[10, 2]], "quad8 reduced");
+ assert!(q8r[0] < 0.005, "{q8r:?}");
+ let im = study(ElementKind::Quad4, incompatible(), &[[10, 1], [20, 2]], "quad4 incompatible");
+ assert!(im[0] < 0.03 && im[1] < 0.01, "{im:?}");
+ let t6 = study(ElementKind::Tri6, full(), &[[10, 2]], "tri6");
+ assert!(t6[0] < 0.01, "{t6:?}");
+ let meshes = [[5, 1], [10, 2], [20, 4], [40, 8]];
+ for (kind, name) in [(ElementKind::Quad4, "quad4"), (ElementKind::Tri3, "tri3")] {
+ let e = study(kind, full(), &meshes, name);
+ let rates: Vec = e.windows(2).map(|w| (w[0] / w[1]).log2()).collect();
+ eprintln!(" {name} convergence rates {rates:.2?}");
+ assert!(e.windows(2).all(|w| w[1] < w[0]), "{name} converges monotonically");
+ assert!(*rates.last().unwrap() > 1.5, "{name} approaches second order: {rates:?}");
+ }
+}
+
+/// Kirsch's plate with a hole of radius a under remote tension, plane stress, quarter model of
+/// half-width 20a: the hoop stress at the top of the hole tends to 3σ₀.
+#[test]
+fn kirsch_plate_with_hole() {
+ let (a, w, s0) = (1.0f64, 20.0f64, 1.0);
+ let ratio: f64 = 40.0;
+ let map = |q: [f64; 2]| {
+ let th = PI / 2.0 * q[1];
+ let outer = if th <= PI / 4.0 { [w, w * th.tan()] } else { [w / th.tan(), w] };
+ let s = (ratio.powf(q[0]) - 1.0) / (ratio - 1.0);
+ let inner = [a * th.cos(), a * th.sin()];
+ [inner[0] + s * (outer[0] - inner[0]), inner[1] + s * (outer[1] - inner[1])]
+ };
+ let mut results = Vec::new();
+ for (kind, formulation, n, name) in [
+ (ElementKind::Quad8, full(), [6, 8], "quad8 coarse"),
+ (ElementKind::Quad8, full(), [12, 16], "quad8 fine"),
+ (ElementKind::Tri6, full(), [12, 16], "tri6"),
+ (ElementKind::Quad4, incompatible(), [24, 32], "quad4 incompatible"),
+ ] {
+ let mut g = grid(n, kind, map);
+ g.define(SectionKind::PlaneStress { thickness: 0.1 }, 1000.0, 0.3, None, formulation, Procedure::Static {});
+ g.symmetry("x0", Symmetry::X, |q| near(q[1], 1.0));
+ g.symmetry("y0", Symmetry::Y, |q| near(q[1], 0.0));
+ g.edge_set("right", |q| near(q[0], 1.0) && q[1] <= 0.5 + 1e-9);
+ g.load("pull", Load::Traction { set: "right".into(), vector: [s0, 0.0, 0.0] });
+ let r = g.solve();
+ let top = g.nodes_where(|q| near(q[0], 0.0) && near(q[1], 1.0));
+ let kt = top.iter().map(|l| r.stress[at(&r, *l)][0]).sum::() / top.len() as f64 / s0;
+ eprintln!("kirsch {name} {n:?}: Kt = {kt:.4} (error vs 3: {:.2}%)", 100.0 * (kt - 3.0).abs() / 3.0);
+ assert!(r.stress.iter().all(|s| s[2] == 0.0), "plane stress has no σzz");
+ results.push(kt);
+ }
+ assert!((results[1] - 3.0).abs() < (results[0] - 3.0).abs(), "{results:?}");
+ assert!((results[1] - 3.0).abs() / 3.0 < 0.02, "{results:?}");
+ assert!((results[2] - 3.0).abs() / 3.0 < 0.03, "{results:?}");
+ assert!((results[3] - 3.0).abs() / 3.0 < 0.04, "{results:?}");
+}
+
+/// Lamé's thick cylinder a ≤ r ≤ b under internal pressure p in plane strain:
+/// σ_r = A − B/r², σ_θ = A + B/r², σ_z = 2νA, u_r = (1 + ν)/E ((1 − 2ν) A r + B/r),
+/// A = p a²/(b² − a²), B = p a² b²/(b² − a²). Solved as a quarter section in plane strain and as
+/// an axisymmetric slice held axially on both faces.
+#[test]
+fn lame_thick_cylinder() {
+ let (a, b, p, young, nu) = (1.0f64, 2.0f64, 10.0, 1000.0, 0.3);
+ let big_a = p * a * a / (b * b - a * a);
+ let big_b = p * a * a * b * b / (b * b - a * a);
+ let ur = |r: f64| (1.0 + nu) / young * ((1.0 - 2.0 * nu) * big_a * r + big_b / r);
+ let sigma_t = |r: f64| big_a + big_b / (r * r);
+ for (kind, formulation, n, name) in [
+ (ElementKind::Quad8, full(), [6, 8], "quad8"),
+ (ElementKind::Quad4, incompatible(), [12, 16], "quad4 incompatible"),
+ (ElementKind::Tri6, full(), [6, 8], "tri6"),
+ ] {
+ // Plane strain quarter: (r, θ) mapped.
+ let mut g = grid(n, kind, |q| {
+ let (r, th) = (a + (b - a) * q[0], PI / 2.0 * q[1]);
+ [r * th.cos(), r * th.sin()]
+ });
+ g.define(STRAIN, young, nu, None, formulation, Procedure::Static {});
+ g.symmetry("y0", Symmetry::Y, |q| near(q[1], 0.0));
+ g.symmetry("x0", Symmetry::X, |q| near(q[1], 1.0));
+ g.edge_set("bore", |q| near(q[0], 0.0));
+ g.load("p", Load::Pressure { set: "bore".into(), magnitude: p });
+ let r = g.solve();
+ let mut eu: f64 = 0.0;
+ for (label, u) in r.nodes.iter().zip(&r.displacement) {
+ let x = g.model.node(*label).unwrap().xyz;
+ let rad = x[0].hypot(x[1]);
+ eu = eu.max(((u[0] * x[0] + u[1] * x[1]) / rad - ur(rad)).abs() / ur(a));
+ }
+ let bore = g.nodes_where(|q| near(q[0], 0.0) && near(q[1], 0.5));
+ let st = bore
+ .iter()
+ .map(|l| {
+ let s = r.stress[at(&r, *l)];
+ let x = g.model.node(*l).unwrap().xyz;
+ let (c, sn) = (x[0] / x[0].hypot(x[1]), x[1] / x[0].hypot(x[1]));
+ s[0] * sn * sn - 2.0 * s[3] * sn * c + s[1] * c * c
+ })
+ .sum::()
+ / bore.len() as f64;
+ let es = (st - sigma_t(a)).abs() / sigma_t(a);
+ eprintln!("lamé plane strain {name} {n:?}: u_r error {:.3}%, bore hoop {st:.4} vs {:.4} ({:.2}%)", 100.0 * eu, sigma_t(a), 100.0 * es);
+ assert!(eu < 0.01 && es < 0.03, "{name}");
+
+ // Axisymmetric slice r ∈ [a, b], 0 ≤ y ≤ 0.5, held axially on both faces.
+ let mut g = grid([n[1], 2], kind, |q| [a + (b - a) * q[0], 0.5 * q[1]]);
+ g.define(AXI, young, nu, None, formulation, Procedure::Static {});
+ g.symmetry("ends", Symmetry::Y, |q| near(q[1], 0.0) || near(q[1], 1.0));
+ g.edge_set("bore", |q| near(q[0], 0.0));
+ g.load("p", Load::Pressure { set: "bore".into(), magnitude: p });
+ let r = g.solve();
+ let mut eu: f64 = 0.0;
+ for (label, u) in r.nodes.iter().zip(&r.displacement) {
+ let x = g.model.node(*label).unwrap().xyz;
+ eu = eu.max((u[0] - ur(x[0])).abs() / ur(a));
+ }
+ let bore = g.nodes_where(|q| near(q[0], 0.0));
+ let st = bore.iter().map(|l| r.stress[at(&r, *l)][2]).sum::() / bore.len() as f64;
+ let sz = r.stress.iter().zip(&r.nodes).map(|(s, _)| s[1]).sum::() / r.nodes.len() as f64;
+ let es = (st - sigma_t(a)).abs() / sigma_t(a);
+ eprintln!(
+ "lamé axisymmetric {name}: u_r error {:.3}%, bore hoop {st:.4} vs {:.4} ({:.2}%), mean σ_z {sz:.4} vs {:.4}",
+ 100.0 * eu,
+ sigma_t(a),
+ 100.0 * es,
+ 2.0 * nu * big_a
+ );
+ assert!(eu < 0.01 && es < 0.03 && (sz - 2.0 * nu * big_a).abs() < 0.02 * sigma_t(a), "{name}");
+ }
+}
+
+/// Rotating discs and cylinders (S. Timoshenko, J. N. Goodier, *Theory of Elasticity*, 3rd ed.,
+/// McGraw-Hill 1970, §32), solid, radius b, speed ω, density ρ:
+/// - a thin disc (plane stress, rotating about Z): σ_r = (3 + ν)/8 ρω² (b² − r²),
+/// σ_θ = ρω²/8 ((3 + ν) b² − (1 + 3ν) r²);
+/// - a long cylinder (axisymmetric, axial strain held at zero):
+/// σ_r = (3 − 2ν)/(8(1 − ν)) ρω² (b² − r²), σ_θ = ρω²/(8(1 − ν)) ((3 − 2ν) b² − (1 + 2ν) r²).
+#[test]
+fn rotating_solid_disc_and_cylinder() {
+ let (b, rho, w, young, nu) = (1.0f64, 1e-3, 30.0, 1000.0, 0.3);
+ let k = rho * w * w;
+ let disc = |r: f64| [(3.0 + nu) / 8.0 * k * (b * b - r * r), k / 8.0 * ((3.0 + nu) * b * b - (1.0 + 3.0 * nu) * r * r)];
+ let cylinder = |r: f64| {
+ let c = k / (8.0 * (1.0 - nu));
+ [c * (3.0 - 2.0 * nu) * (b * b - r * r), c * ((3.0 - 2.0 * nu) * b * b - (1.0 + 2.0 * nu) * r * r)]
+ };
+ let peak = disc(0.0)[0];
+ // Linear elements resolve the quadratic stress field to about their size (h/b of the peak
+ // in nodal averages); they get a finer mesh.
+ for (kind, formulation, n, tol, name) in [
+ (ElementKind::Quad8, full(), 8, 0.01, "quad8"),
+ (ElementKind::Tri6, full(), 8, 0.01, "tri6"),
+ (ElementKind::Quad4, incompatible(), 32, 0.04, "quad4 incompatible"),
+ ] {
+ // Plane stress quarter disc, polar-mapped from a pin hole of radius b/50 (a polar map
+ // folds at the centre); the hole's disturbance decays as (r₀/r)², so the solid disc's
+ // field is compared at r ≥ b/2.
+ let mut g = grid([n, 8], kind, |q| {
+ let (r, th) = (b * (0.02 + 0.98 * q[0]), PI / 2.0 * q[1]);
+ [r * th.cos(), r * th.sin()]
+ });
+ g.define(SectionKind::PlaneStress { thickness: 0.05 }, young, nu, Some(rho), formulation, Procedure::Static {});
+ g.symmetry("y0", Symmetry::Y, |q| near(q[1], 0.0));
+ g.symmetry("x0", Symmetry::X, |q| near(q[1], 1.0));
+ g.load("spin", Load::Centrifugal { angular_velocity: w, point: [0.0; 3], axis: [0.0, 0.0, 1.0], set: None });
+ let r = g.solve();
+ let mut e: f64 = 0.0;
+ for (label, s) in r.nodes.iter().zip(&r.stress) {
+ let x = g.model.node(*label).unwrap().xyz;
+ let rad = x[0].hypot(x[1]);
+ if rad < 0.5 * b {
+ continue;
+ }
+ let (c, sn) = (x[0] / rad, x[1] / rad);
+ let sr = s[0] * c * c + 2.0 * s[3] * sn * c + s[1] * sn * sn;
+ let st = s[0] * sn * sn - 2.0 * s[3] * sn * c + s[1] * c * c;
+ let want = disc(rad);
+ e = e.max((sr - want[0]).abs().max((st - want[1]).abs()) / peak);
+ }
+ eprintln!("rotating disc (plane stress) {name}: stress error {:.3}% of the peak (r ≥ b/2)", 100.0 * e);
+ assert!(e < tol, "{name}: {e}");
+
+ // Axisymmetric solid cylinder 0 ≤ r ≤ b, a slice held axially on both faces: nodes on the
+ // axis are held radially by the solver.
+ let mut g = grid([n, 2], kind, |q| [b * q[0], 0.2 * q[1]]);
+ g.define(AXI, young, nu, Some(rho), formulation, Procedure::Static {});
+ g.symmetry("ends", Symmetry::Y, |q| near(q[1], 0.0) || near(q[1], 1.0));
+ g.load("spin", Load::Centrifugal { angular_velocity: w, point: [0.0; 3], axis: [0.0, 1.0, 0.0], set: None });
+ let r = g.solve();
+ let cyl_peak = cylinder(0.0)[0];
+ let mut e: f64 = 0.0;
+ for (label, s) in r.nodes.iter().zip(&r.stress) {
+ let x = g.model.node(*label).unwrap().xyz[0];
+ let want = cylinder(x);
+ e = e.max((s[0] - want[0]).abs().max((s[2] - want[1]).abs()) / cyl_peak);
+ }
+ let on_axis = g.nodes_where(|q| near(q[0], 0.0));
+ assert!(on_axis.iter().all(|l| r.displacement[at(&r, *l)][0] == 0.0), "the axis does not move radially");
+ eprintln!("rotating cylinder (axisymmetric) {name}: stress error {:.3}% of the peak", 100.0 * e);
+ assert!(e < 2.0 * tol, "{name}: {e}");
+ }
+}
+
+/// The total mass of a 2D model is its area times its thickness, or the volume of revolution
+/// (Pappus) when axisymmetric; a free plate has three rigid-body modes in its plane.
+#[test]
+fn masses_and_rigid_modes() {
+ let rho = 2.0;
+ for (theory, mass) in [(STRESS, rho * 2.0 * 0.5), (AXI, rho * 2.0 * PI * 1.5 * 2.0)] {
+ let mut g = grid([4, 2], ElementKind::Quad8, |p| [1.0 + p[0], 2.0 * p[1]]);
+ g.define(theory, 1000.0, 0.3, Some(rho), full(), Procedure::Frequency { modes: 4, min_hz: None, max_hz: None });
+ let f = solve_frequency(&g.model, None, &FrequencyOptions::default(), &Limits::default()).unwrap();
+ eprintln!("{}: mass {} vs {mass}, frequencies {:?}", theory.name(), f.total_mass, f.modes.iter().map(|m| m.frequency).collect::>());
+ assert!((f.total_mass - mass).abs() < 1e-12 * mass);
+ let rigid = f.modes.iter().filter(|m| m.frequency < 1e-3 * f.modes.last().unwrap().frequency).count();
+ // In-plane: two translations and a rotation; axisymmetric: the axial translation only.
+ assert_eq!(rigid, if theory == AXI { 1 } else { 3 }, "{}", theory.name());
+ }
+}
+
+#[test]
+fn two_dimensional_results_are_deterministic_and_errors_name_the_problem() {
+ let make = || {
+ let mut g = grid([4, 2], ElementKind::Quad4, |p| [1.0 + p[0] + 0.1 * p[1] * p[0], p[1]]);
+ g.define(AXI, 1000.0, 0.3, Some(1.0), incompatible(), Procedure::Static {});
+ g.symmetry("bottom", Symmetry::Y, |p| near(p[1], 0.0));
+ g.edge_set("bore", |p| near(p[0], 0.0));
+ g.load("p", Load::Pressure { set: "bore".into(), magnitude: 3.0 });
+ g.load("g", Load::Gravity { acceleration: [0.0, -9.81, 0.0], set: None });
+ g
+ };
+ let g = make();
+ assert_eq!(g.solve(), g.solve());
+ // A coupling is not supported in 2D yet.
+ let mut g = make();
+ g.model.add_reference_point("RP", modelcraft_model::ReferencePoint { xyz: [0.0, 2.0, 0.0], label: 0 }).unwrap();
+ g.model
+ .add_coupling(
+ "c",
+ modelcraft_model::Coupling {
+ reference_point: "RP".into(),
+ set: "bottom".into(),
+ kind: modelcraft_model::CouplingKind::Kinematic,
+ dofs: modelcraft_model::Dofs::all(),
+ },
+ )
+ .unwrap();
+ let e = solve_static(&g.model, None, &Limits::default()).unwrap_err().to_string();
+ assert!(e.contains("couplings in 2D models"), "{e}");
+ // A load with a Z component, and a rotation about a wrong axis, are refused by the check.
+ let mut g = make();
+ g.load("lift", Load::Gravity { acceleration: [0.0, 0.0, 1.0], set: None });
+ let e = solve_static(&g.model, None, &Limits::default()).unwrap_err().to_string();
+ assert!(e.contains("Z component"), "{e}");
+ // No support: singular, naming a node and direction.
+ let mut g = grid([2, 1], ElementKind::Quad4, |p| [p[0], p[1]]);
+ g.define(STRESS, 1000.0, 0.3, None, full(), Procedure::Static {});
+ g.edge_set("right", |p| near(p[0], 1.0));
+ g.load("p", Load::Pressure { set: "right".into(), magnitude: 1.0 });
+ let e = solve_static(&g.model, None, &Limits::default()).unwrap_err().to_string();
+ assert!(e.contains("singular") && e.contains("node"), "{e}");
+}
diff --git a/docs/architecture.md b/docs/architecture.md
index fe7a894..923f025 100644
--- a/docs/architecture.md
+++ b/docs/architecture.md
@@ -133,14 +133,19 @@ module docs; ADR 0001 records why the linear algebra is our own.
1. **System.** `Model::resolve()`; each solid section's elements take their material's
elasticity matrix and the section's formulation (`modelcraft_fem::element::ELEMENT_TYPES`, one
- row per element kind: tet4, tet10, wedge6, hex8, hex20). Every reference point a coupling
+ row per element kind: tet4, tet10, wedge6, hex8, hex20); 2D sections (plane stress, plane
+ strain, axisymmetric) take the 2D table (`modelcraft_fem::plane::PLANE_TYPES`: tri3, tri6,
+ quad4, quad8), chosen by the section kind, and the model then has two degrees of freedom per
+ node: the 2D elements declare `ElementDofs::new(DofSet::PLANE)` (Ux, Uy) to the `DofLayout`. A 2D model
+ lies in the XY plane (axisymmetric: x the radius, Y the axis; nodes on the axis are held
+ radially), is loaded on element edges, and cannot mix with solids or couplings yet. Every reference point a coupling
ties to the mesh gets six degrees of freedom (translations, rotations) after the mesh's. Shells
and beams return "not supported yet". The step's boundary conditions prescribe translations
(rotations do not exist on solid nodes) and a coupled point's six directions; its loads become
a load vector: forces and tractions from `Resolved::nodal_loads`, forces and moments on coupled
reference points, pressure as consistent face loads, gravity and body forces through the
element shape functions.
-2. **Assembly** on the node graph (3 × 3 blocks, elements in label order, so sums are
+2. **Assembly** on the node graph (3 × 3 blocks for solids, 2 × 2 in 2D models, elements in label order, so sums are
reproducible bit for bit).
3. **Ordering and elimination** (`constraints`). Nested dissection of the nodes by coordinate
bisection. Couplings become linear constraint rows: kinematic = rigid link of each coupled
@@ -167,6 +172,11 @@ Verification (`crates/solver/tests/verification.rs`, analytic only): displacemen
patch tests for every element and formulation (exact to round-off), cantilever tip deflection
against Timoshenko with convergence rates of the linear elements, Lamé's cylinder, Kirsch's
plate with a hole (Kt → 3), cantilever and free–free beam frequencies with six rigid-body modes.
+`crates/solver/tests/plane.rs`: the 2D elements' displacement patch tests in plane stress, plane
+strain and axisymmetry, uniform traction, a plane stress cantilever against Timoshenko with
+convergence rates, Kirsch's plate, Lamé's cylinder in plane strain and axisymmetric, a rotating
+disc (plane stress) and a rotating cylinder (axisymmetric) against Timoshenko & Goodier, masses
+and rigid-body modes.
`crates/solver/tests/couplings.rs`: a block given a rigid motion through a coupled point
(kinematic and distributing), clamping through a coupled point equals clamping (static and
modes), tension through a distributing coupling against the equivalent pressure away from the
diff --git a/docs/features.md b/docs/features.md
index f4b506c..eb2e996 100644
--- a/docs/features.md
+++ b/docs/features.md
@@ -5,7 +5,7 @@ Generated from the capability catalogue (`crates/engine/src/catalogue.rs`) by
Every entry is a toolbar button. **Ready** entries run their command; **coming soon** entries are
shown greyed with the roadmap milestone that brings them (ROADMAP.md).
-45 ready, 55 coming soon.
+46 ready, 55 coming soon.
## Geometry
@@ -30,6 +30,7 @@ shown greyed with the roadmap milestone that brings them (ROADMAP.md).
| Tet | Tet Mesh — Constrained Delaunay tet mesh of a body: every boundary triangle kept, faces tagged (refinement comes next) | ready | `mesh.tet.generate` |
| Tet | Quadratic — Second-order elements with mid-side nodes on the geometry | coming soon (M2) | |
| Tet | Delaunay Points — Delaunay tetrahedralization of a point list (debugging aid) | ready | `mesh.delaunay.points` |
+| 2D | 2D Mesh — Quad or triangle mesh of a planar face or polygon in the XY plane for plane stress, plane strain and axisymmetric models; edges tagged | ready | `mesh.planar` |
| Structured Hex | Quad Faces — All-quad meshes of faces on shared edge nodes | coming soon (M5) | |
| Structured Hex | Mapped Hex — Mapped hex mesh of six-sided regions | coming soon (M5) | |
| Structured Hex | Sweep Hex — Hex mesh swept along extrusions and revolutions | coming soon (M5) | |
@@ -60,7 +61,7 @@ shown greyed with the roadmap milestone that brings them (ROADMAP.md).
| Materials | Damage — Damage initiation and evolution | coming soon (M8) | |
| Materials | Thermal — Conductivity, specific heat and expansion | coming soon (M8) | |
| Materials | Equation of State — Polynomial, ideal gas, Mie–Grüneisen, JWL and tabulated equations of state | coming soon (M11) | |
-| Sections | Section — Solid, shell, beam or rod section assigned to a set, with its element formulation | ready | `model.section.create`, `model.section.edit`, `model.section.delete` |
+| Sections | Section — Solid, shell, beam, rod or 2D (plane stress, plane strain, axisymmetric) section assigned to a set, with its element formulation | ready | `model.section.create`, `model.section.edit`, `model.section.delete` |
| Sections | Composite — Composite layups and continuum shells | coming soon (M8) | |
| Sections | Orientation — Material orientations | coming soon (M8) | |
| Sets | Set — A set of nodes, elements and faces, or of named geometry that survives remeshing | ready | `model.set.create`, `model.set.edit`, `model.set.delete`, `model.set.resolve` |
diff --git a/docs/meshing/literature.md b/docs/meshing/literature.md
index f428ced..68db8c9 100644
--- a/docs/meshing/literature.md
+++ b/docs/meshing/literature.md
@@ -849,6 +849,13 @@ boundary nodes, so faces sharing an edge share its nodes. DOIs below were checke
[10.1016/S0925-7721(01)00047-5](https://doi.org/10.1016/S0925-7721(01)00047-5). Used without segment splitting
(boundary nodes are prescribed).
+**B7.3a Free triangle meshes on prescribed boundary nodes** (`free_tri_mesh`, used by
+`mesh.planar` for tri3/tri6 2D models): the same constrained Delaunay triangulation of the full
+resolution loops, refined to the size by circumcentre insertion that never splits a boundary
+segment (J. Ruppert, *A Delaunay refinement algorithm for quality 2-dimensional mesh generation*,
+J. Algorithms 18(3):548–585, 1995, without segment splitting, so the boundary nodes stay the
+input points), Steiner points kept half the local spacing from the boundary.
+
**B7.4 Smoothing and untangling of the quads.** **Verdict: BUILD ON.**
- L. A. Freitag. *On combining Laplacian and optimization-based mesh smoothing techniques.* AMD-Vol. 220, ASME 1997;
S. A. Canann, J. R. Tristano, M. L. Staten, *An approach to combined Laplacian and optimization-based smoothing for
diff --git a/docs/parity.md b/docs/parity.md
index 979fc82..ae545bc 100644
--- a/docs/parity.md
+++ b/docs/parity.md
@@ -33,7 +33,7 @@ products are run, no documentation text is copied). Status: **done**, **started*
| Mesh | Size field: global, per face and edge seeds, bias, curvature and proximity refinement | M2 | yes | planned |
| Mesh | Free tet meshing (tet4, tet10 with mid-side nodes on the geometry) | M2 | | started |
| Mesh | Quality metrics, certified validity, verify per set | M1–M2 | yes | done (metrics, gate) |
-| Mesh | 2D quad meshing of faces | M5 | yes | planned |
+| Mesh | 2D quad meshing of faces | M5 | yes | started (`modelcraft_mesh::quad`; `mesh.planar`: 2D models of a polygon or a planar face, quad4/quad8/tri3/tri6, edges tagged by their neighbouring face) |
| Mesh | Mapped hex of six-sided regions | M5 | yes | planned |
| Mesh | Sweep hex (extrusions, revolutions) | M5 | yes | planned |
| Mesh | Automatic decomposition into sweepable and mappable regions | M5 | yes | planned |
@@ -46,6 +46,7 @@ products are run, no documentation text is copied). Status: **done**, **started*
| 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 | | done |
+| Model | 2D sections: plane stress, plane strain, axisymmetric (XY plane, Y the axis); element edges in sets; centrifugal load | M4 | | 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 | done (model; per-solver mapping with the deck writers) |
| Model | Boundary conditions: fixed, displacement, symmetry | M3 | | done |
@@ -60,7 +61,7 @@ products are run, no documentation text is copied). Status: **done**, **started*
| Decks | Nastran bulk data writer (`deck.write`, SOL 101 or 103 per deck; reports what Nastran cannot express) | M3 | | done |
| Decks | LS-DYNA keyword writer (implicit; explicit in M8) | M3 | | done (implicit static and eigenvalue; beam sections planned) |
| Decks | Readers keeping unknown keywords | M9 | | planned |
-| Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | started (solid elements: tet4, tet10, wedge6, hex8, hex20; reference points, kinematic and distributing couplings, moments) |
+| Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | started (solid elements: tet4, tet10, wedge6, hex8, hex20; 2D plane stress, plane strain and axisymmetric tri3, tri6, quad4 (full, incompatible modes), quad8 (full, reduced); reference points, kinematic and distributing couplings, moments; centrifugal loads) |
| Solve | Own explicit dynamics solver (central differences, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, energy balance), verified against analytic wave and impact solutions | M10 | yes | started (X1 core: small-strain elastic, rigid walls; large deformation, plasticity and contact planned) |
| Results | Deformed shape, contours, mode shapes, probes, reaction sums | M4 | | started (results in the session, `solve.probe`, reaction sums; display planned) |
| Results | XY and path plots, cuts, display groups, free-body forces, animation | M9 | | planned |
@@ -98,6 +99,15 @@ products are run, no documentation text is copied). Status: **done**, **started*
body forces per volume, partial-DOF kinematic couplings, constraints on a coupled reference
point, gravity on different sets or on parts of sections, hybrid other than tet4, reduced
quadratic elements, materials without density, more than 10 shell integration points.
+- 2D models: Abaqus writes CPS*/CPE*/CAX* (3, 6, 4, 4R, 4I, 8, 8R) with the thickness on the
+ `*SOLID SECTION` line, edge surfaces `S1`… and `BR`/`BZ` body forces for axisymmetric elements;
+ Nastran writes plane stress as membranes (`PSHELL` with MID1 only, grids holding 3456 through
+ `GRID` PS) and edge pressures as consistent `FORCE`s, and refuses plane strain and
+ axisymmetric elements (their cards differ between Nastran dialects); LS-DYNA writes plane stress
+ and plane strain (unit thickness) as `*SECTION_SHELL` ELFORM 12/13, one-point quad4 (reduced
+ integration, with hourglass control) and tri3, and refuses quadratic 2D elements, full or
+ incompatible-mode quad4 and axisymmetric sections. Centrifugal loads: Abaqus `CENTRIF`
+ (ω², point, axis), Nastran `RFORCE` about an axis through the origin, LS-DYNA not yet.
- Open item (LS-DYNA tet10 node order): the writer follows the manual's figure (mid-side nodes 8,
9, 10 on edges 2-4, 3-4, 1-4), which disagrees with an earlier reading as the Abaqus order. It
needs one confirmation run by a licensed user; agents never run commercial solvers.
diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md
index 80fa106..b0c305f 100644
--- a/docs/solvers/literature.md
+++ b/docs/solvers/literature.md
@@ -979,6 +979,14 @@ the module docs):
Lewis, H. D. Simon, SIAM J. Matrix Anal. Appl. 15(1):228–272, 1994; B. N. Parlett, *The
Symmetric Eigenvalue Problem*, SIAM 1998; G. H. Golub, C. F. Van Loan, *Matrix Computations*,
4th ed., 2013.
+- 2D elements (plane stress, plane strain, axisymmetric): Zienkiewicz, Taylor & Zhu 2013, ch. 6
+ and §6.6; Hughes 2000, §2.7 and §3.5. Quad4 incompatible modes (QM6): Wilson et al. 1973 and
+ Taylor, Beresford & Wilson 1976 as for hex8; in axisymmetric elements the modes' strain
+ operator is shifted by its volume mean, the patch-test correction of E. L. Wilson and A.
+ Ibrahimbegovic, *Use of incompatible displacement modes for the calculation of element
+ stiffnesses or stresses*, Finite Elements in Analysis and Design 7(3):229–241, 1990. Rotating
+ discs and cylinders for verification: S. Timoshenko, J. N. Goodier, *Theory of Elasticity*,
+ 3rd ed., McGraw-Hill 1970, §32.
- Wedge6: the linear triangle times the linear interval; the symmetric three-point triangle rule
(Zienkiewicz, Taylor & Zhu 2013, table 6.3) times two-point Gauss–Legendre.
- Couplings and constraint elimination: rigid links and master–slave elimination by