diff --git a/crates/specialized/rtx-fea/src/analysis/flag3d.rs b/crates/specialized/rtx-fea/src/analysis/flag3d.rs new file mode 100644 index 0000000..9af8e57 --- /dev/null +++ b/crates/specialized/rtx-fea/src/analysis/flag3d.rs @@ -0,0 +1,475 @@ +//! The Turek–Hron flag as a 3-D solid: a Hex20 plate `[x0, x1] × [y0, y1] +//! × [z0, z1]` (length × thickness × span), root face `x = x0` clamped, +//! for R8's coupled 3-D FSI (omni-cortex roadmap, item R8-b). +//! +//! Nothing here is a new solver: [`NonlinearDynamicAnalysis`] and the +//! total-Lagrangian St. Venant–Kirchhoff path are dimension-generic, so +//! the 3-D flag steps through exactly the machinery the 2-D harness uses +//! (`set_nodal_forces`, `step(state) → state`, Newmark γ/β, the +//! line-search/subdivision rescue). This module adds what the 2-D harness +//! builds by hand: +//! +//! * the structured Hex20 mesh on the `(2nx+1) × (2ny+1) × (2nz+1)` +//! serendipity lattice (x-major node order, so the sequential DOF +//! numbering keeps the banded LU's bandwidth at one x-slab); +//! * the root clamp, with the lateral faces either **free** (the real +//! plate, the physical 3-D problem) or **plane strain** (`u_z = 0` at +//! every node: a z-independent field is exactly representable by Hex20 +//! and the 3-D energy then equals span × the 2-D plane-strain energy, so +//! this reproduces the 2-D Quad8 model to rounding — the self-consistency +//! gate, "run the 3-D problem as the 2-D problem first"); +//! * the wetted surface (bottom, top, tip, and the two lateral faces) as +//! Quad8 faces with outward orientation, and the consistent nodal forces +//! of a traction field integrated over the *current* (deformed) faces — +//! the interface load a partitioned coupling feeds to +//! [`NonlinearDynamicStepper::set_nodal_forces`](super::NonlinearDynamicStepper::set_nodal_forces). +//! +//! New code only: no existing solver path changes. + +use super::{AnalysisConfig, ConvergenceCriteria, NonlinearDynamicAnalysis}; +use crate::assembly::dof_mapping::DofComponent; +use crate::boundary::dirichlet::{DirichletBC, DirichletType}; +use crate::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction}; +use crate::error::FeaResult; +use crate::materials::{LinearElastic, MaterialDatabase}; +use crate::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId}; +use nalgebra::{DVector, Vector3}; +use std::collections::BTreeMap; + +/// Geometry and resolution of the plate. `nx`, `ny`, `nz` are Hex20 +/// element counts along length, thickness and span. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Flag3dSpec { + pub x0: f64, + pub x1: f64, + pub y0: f64, + pub y1: f64, + pub z0: f64, + pub z1: f64, + pub nx: usize, + pub ny: usize, + pub nz: usize, +} + +impl Flag3dSpec { + /// The Turek–Hron flag (`[0.25, 0.6] × [0.19, 0.21]`) extruded over + /// `z ∈ [z0, z0 + span]`. + pub fn turek_hron(span: f64, z0: f64, nx: usize, ny: usize, nz: usize) -> Self { + Self { + x0: 0.25, + x1: 0.6, + y0: 0.19, + y1: 0.21, + z0, + z1: z0 + span, + nx, + ny, + nz, + } + } + + pub fn span(&self) -> f64 { + self.z1 - self.z0 + } +} + +/// How the lateral faces `z = z0, z1` are held. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum LateralFaces { + /// Traction-free: the physical 3-D plate. + Free, + /// `u_z = 0` at every node: exact 2-D plane strain (the 2-D model's + /// definition), for the self-consistency gate. + PlaneStrain, +} + +/// Which wetted face of the plate a [`SurfaceFace`] lies on. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub enum FlagSide { + /// `y = y0` + Bottom, + /// `y = y1` + Top, + /// `x = x1` + Tip, + /// `z = z0` + SideLow, + /// `z = z1` + SideHigh, +} + +/// A Quad8 face of the Hex20 mesh: corners counter-clockwise seen from +/// outside, then mid-edge nodes `(0-1, 1-2, 2-3, 3-0)` — the library's +/// Quad8 order, oriented so `∂x/∂ξ × ∂x/∂η` points out of the solid. +#[derive(Debug, Clone, PartialEq)] +pub struct SurfaceFace { + pub side: FlagSide, + pub nodes: [NodeId; 8], +} + +/// The structured Hex20 plate. +#[derive(Debug, Clone)] +pub struct Flag3d { + pub spec: Flag3dSpec, + pub mesh: Mesh, + lattice: Vec>, + dims: [usize; 3], +} + +const GAUSS3: [(f64, f64); 3] = [ + (-0.774_596_669_241_483_4, 5.0 / 9.0), + (0.0, 8.0 / 9.0), + (0.774_596_669_241_483_4, 5.0 / 9.0), +]; + +/// Quad8 serendipity shape functions and their `(ξ, η)` derivatives, the +/// library's node order. +fn quad8(xi: f64, eta: f64) -> ([f64; 8], [[f64; 2]; 8]) { + let corners = [(-1.0, -1.0), (1.0, -1.0), (1.0, 1.0), (-1.0, 1.0)]; + let mut n = [0.0; 8]; + let mut d = [[0.0; 2]; 8]; + for (a, &(xa, ya)) in corners.iter().enumerate() { + let (p, q) = (1.0 + xa * xi, 1.0 + ya * eta); + let r = xa * xi + ya * eta - 1.0; + n[a] = 0.25 * p * q * r; + d[a][0] = 0.25 * xa * (q * r + p * q); + d[a][1] = 0.25 * ya * (p * r + p * q); + } + // Mid-edge nodes: (0, -1), (1, 0), (0, 1), (-1, 0). + n[4] = 0.5 * (1.0 - xi * xi) * (1.0 - eta); + d[4] = [-xi * (1.0 - eta), -0.5 * (1.0 - xi * xi)]; + n[5] = 0.5 * (1.0 + xi) * (1.0 - eta * eta); + d[5] = [0.5 * (1.0 - eta * eta), -(1.0 + xi) * eta]; + n[6] = 0.5 * (1.0 - xi * xi) * (1.0 + eta); + d[6] = [-xi * (1.0 + eta), 0.5 * (1.0 - xi * xi)]; + n[7] = 0.5 * (1.0 - xi) * (1.0 - eta * eta); + d[7] = [-0.5 * (1.0 - eta * eta), -(1.0 - xi) * eta]; + (n, d) +} + +impl Flag3d { + /// Build the mesh (material 0 on every element). + pub fn build(spec: Flag3dSpec) -> FeaResult { + assert!( + spec.nx > 0 && spec.ny > 0 && spec.nz > 0, + "element counts must be positive" + ); + let dims = [2 * spec.nx + 1, 2 * spec.ny + 1, 2 * spec.nz + 1]; + let mut mesh = Mesh::new(3)?; + let mut lattice = vec![None; dims[0] * dims[1] * dims[2]]; + for i in 0..dims[0] { + for j in 0..dims[1] { + for k in 0..dims[2] { + if (i % 2) + (j % 2) + (k % 2) > 1 { + continue; // not a serendipity node + } + let x = spec.x0 + (spec.x1 - spec.x0) * i as f64 / (dims[0] - 1) as f64; + let y = spec.y0 + (spec.y1 - spec.y0) * j as f64 / (dims[1] - 1) as f64; + let z = spec.z0 + (spec.z1 - spec.z0) * k as f64 / (dims[2] - 1) as f64; + lattice[(i * dims[1] + j) * dims[2] + k] = + Some(mesh.add_node(Node::new_3d(x, y, z))); + } + } + } + let mut flag = Self { + spec, + mesh, + lattice, + dims, + }; + for ex in 0..spec.nx { + for ey in 0..spec.ny { + for ez in 0..spec.nz { + let (a, b, c) = (2 * ex, 2 * ey, 2 * ez); + let at = |i, j, k| flag.lattice_node(i, j, k).expect("serendipity node"); + let nodes = vec![ + at(a, b, c), + at(a + 2, b, c), + at(a + 2, b + 2, c), + at(a, b + 2, c), + at(a, b, c + 2), + at(a + 2, b, c + 2), + at(a + 2, b + 2, c + 2), + at(a, b + 2, c + 2), + at(a + 1, b, c), + at(a + 2, b + 1, c), + at(a + 1, b + 2, c), + at(a, b + 1, c), + at(a + 1, b, c + 2), + at(a + 2, b + 1, c + 2), + at(a + 1, b + 2, c + 2), + at(a, b + 1, c + 2), + at(a, b, c + 1), + at(a + 2, b, c + 1), + at(a + 2, b + 2, c + 1), + at(a, b + 2, c + 1), + ]; + flag.mesh.add_element(Element::new( + ElementType::Hex20, + nodes, + MaterialId(0), + )?)?; + } + } + } + Ok(flag) + } + + /// The node at lattice index `(i, j, k)` (`0..=2n` per direction), if + /// that lattice point carries a serendipity node. + pub fn lattice_node(&self, i: usize, j: usize, k: usize) -> Option { + if i >= self.dims[0] || j >= self.dims[1] || k >= self.dims[2] { + return None; + } + self.lattice[(i * self.dims[1] + j) * self.dims[2] + k] + } + + /// Lattice sizes `(2nx+1, 2ny+1, 2nz+1)`. + pub fn lattice_dims(&self) -> [usize; 3] { + self.dims + } + + /// The node nearest to `p` (reference coordinates). + pub fn nearest_node(&self, p: Vector3) -> NodeId { + self.mesh + .nodes + .iter() + .min_by(|a, b| { + let da = (a.1.position() - p).norm(); + let db = (b.1.position() - p).norm(); + da.partial_cmp(&db).unwrap() + }) + .map(|(&id, _)| id) + .expect("non-empty mesh") + } + + /// The Turek–Hron point A `(x1, (y0+y1)/2)` on the mid-span line + /// (the nearest lattice node: exact for even `nz`, else the nearer of + /// the two mid-span candidates). + pub fn point_a(&self) -> NodeId { + let s = &self.spec; + self.nearest_node(Vector3::new(s.x1, 0.5 * (s.y0 + s.y1), 0.5 * (s.z0 + s.z1))) + } + + /// Nodes on the clamped root face `x = x0`. + pub fn root_nodes(&self) -> Vec { + let mut out = Vec::new(); + for j in 0..self.dims[1] { + for k in 0..self.dims[2] { + if let Some(id) = self.lattice_node(0, j, k) { + out.push(id); + } + } + } + out + } + + /// The root clamp (`u = 0` on `x = x0`) plus the lateral condition. + pub fn clamp_root(&self, lateral: LateralFaces) -> BoundaryConditionSet { + let zero = || DirichletType::Spatial(SpatialFunction(Box::new(|_| 0.0))); + let dirichlet = |nodes: Vec, component| { + BoundaryCondition::Dirichlet(DirichletBC { + nodes, + components: vec![component], + condition_type: zero(), + time_range: None, + ramping_factor: 1.0, + gradual_enforcement: false, + }) + }; + let root = self.root_nodes(); + let mut set = BoundaryConditionSet::new(); + for component in [ + DofComponent::DisplacementX, + DofComponent::DisplacementY, + DofComponent::DisplacementZ, + ] { + set.add_condition(dirichlet(root.clone(), component)); + } + if lateral == LateralFaces::PlaneStrain { + let all: Vec = self.mesh.nodes.keys().copied().collect(); + set.add_condition(dirichlet(all, DofComponent::DisplacementZ)); + } + set + } + + /// One linear-elastic material (the TL path reads its Lamé pair and + /// density). + pub fn materials(e: f64, nu: f64, rho: f64) -> MaterialDatabase { + let mut db = MaterialDatabase::new(); + db.add_material( + MaterialId(0), + LinearElastic::new(e, nu).with_density(rho), + None, + ); + db + } + + /// The coupled march's structure: total-Lagrangian SVK, Newmark with + /// the given `γ` and `β = (γ + ½)²/4`, 60 Newton iterations — the 2-D + /// harness's settings. `num_steps` only matters for + /// [`NonlinearDynamicAnalysis::run`]; a coupling uses + /// [`NonlinearDynamicAnalysis::stepper`]. + pub fn dynamic_analysis( + &self, + e: f64, + nu: f64, + rho: f64, + lateral: LateralFaces, + dt: f64, + num_steps: usize, + gamma: f64, + ) -> NonlinearDynamicAnalysis { + let beta = (gamma + 0.5).powi(2) / 4.0; + NonlinearDynamicAnalysis::new( + self.mesh.clone(), + Self::materials(e, nu, rho), + self.clamp_root(lateral), + dt, + num_steps, + AnalysisConfig::default(), + ) + .with_total_lagrangian() + .with_convergence_criteria(ConvergenceCriteria { + max_iterations: 60, + ..ConvergenceCriteria::default() + }) + .with_newmark_parameters(gamma, beta) + } + + /// The Quad8 faces of the given sides (the root face is never + /// wetted). Order: side, then element index. + pub fn surface_faces(&self, sides: &[FlagSide]) -> Vec { + let [di, dj, dk] = self.dims; + let (ilast, jlast, klast) = (di - 1, dj - 1, dk - 1); + let mut faces = Vec::new(); + // base lattice point + the two in-face axes (ξ, η) as unit steps. + let mut push = |side, base: [usize; 3], u: [usize; 3], v: [usize; 3]| { + let p = |cu: usize, cv: usize| { + let q = [ + base[0] + cu * u[0] + cv * v[0], + base[1] + cu * u[1] + cv * v[1], + base[2] + cu * u[2] + cv * v[2], + ]; + self.lattice_node(q[0], q[1], q[2]).expect("face node") + }; + faces.push(SurfaceFace { + side, + nodes: [ + p(0, 0), + p(2, 0), + p(2, 2), + p(0, 2), + p(1, 0), + p(2, 1), + p(1, 2), + p(0, 1), + ], + }); + }; + const X: [usize; 3] = [1, 0, 0]; + const Y: [usize; 3] = [0, 1, 0]; + const Z: [usize; 3] = [0, 0, 1]; + for &side in sides { + match side { + FlagSide::Bottom => { + for a in (0..ilast).step_by(2) { + for c in (0..klast).step_by(2) { + push(side, [a, 0, c], X, Z); // x × z = −y + } + } + } + FlagSide::Top => { + for a in (0..ilast).step_by(2) { + for c in (0..klast).step_by(2) { + push(side, [a, jlast, c], Z, X); // z × x = +y + } + } + } + FlagSide::Tip => { + for b in (0..jlast).step_by(2) { + for c in (0..klast).step_by(2) { + push(side, [ilast, b, c], Y, Z); // y × z = +x + } + } + } + FlagSide::SideLow => { + for a in (0..ilast).step_by(2) { + for b in (0..jlast).step_by(2) { + push(side, [a, b, 0], Y, X); // y × x = −z + } + } + } + FlagSide::SideHigh => { + for a in (0..ilast).step_by(2) { + for b in (0..jlast).step_by(2) { + push(side, [a, b, klast], X, Y); // x × y = +z + } + } + } + } + } + faces + } + + /// Every wetted face: bottom, top, tip, and both lateral faces. + pub fn wetted_faces(&self) -> Vec { + self.surface_faces(&[ + FlagSide::Bottom, + FlagSide::Top, + FlagSide::Tip, + FlagSide::SideLow, + FlagSide::SideHigh, + ]) + } + + /// Consistent nodal forces `f_a = ∫ N_a t(x, n) da` of a traction + /// field over the given faces in the configuration `X + u` + /// (`displacement` in the stepper's global DOF numbering, read through + /// `node_dofs`; `None` = the reference configuration). `traction` + /// receives the current point and the current outward unit normal — + /// a pressure `p` is `|x, n| -p(x) * n`. 3 × 3 Gauss per face. + /// Returns one entry per touched node, sorted by node id. + pub fn face_nodal_forces( + &self, + faces: &[SurfaceFace], + displacement: Option<(&DVector, &dyn Fn(NodeId) -> Vec)>, + traction: &dyn Fn(Vector3, Vector3) -> Vector3, + ) -> Vec<(NodeId, Vector3)> { + let position = |id: NodeId| -> Vector3 { + let x = self.mesh.get_node(id).expect("face node").position(); + match displacement { + Some((u, dofs)) => { + let d = dofs(id); + x + Vector3::new(u[d[0]], u[d[1]], u[d[2]]) + } + None => x, + } + }; + let mut out: BTreeMap> = BTreeMap::new(); + for face in faces { + let xs: Vec> = face.nodes.iter().map(|&id| position(id)).collect(); + for &(xi, wx) in &GAUSS3 { + for &(eta, wy) in &GAUSS3 { + let (n, d) = quad8(xi, eta); + let mut x = Vector3::zeros(); + let mut t1 = Vector3::zeros(); + let mut t2 = Vector3::zeros(); + for a in 0..8 { + x += xs[a] * n[a]; + t1 += xs[a] * d[a][0]; + t2 += xs[a] * d[a][1]; + } + let cross = t1.cross(&t2); + let jac = cross.norm(); + let t = traction(x, cross / jac); + let w = wx * wy * jac; + for a in 0..8 { + *out.entry(face.nodes[a]).or_insert_with(Vector3::zeros) += t * (n[a] * w); + } + } + } + } + out.into_iter().collect() + } +} diff --git a/crates/specialized/rtx-fea/src/analysis/mod.rs b/crates/specialized/rtx-fea/src/analysis/mod.rs index 0d33aa4..1dcc2b7 100644 --- a/crates/specialized/rtx-fea/src/analysis/mod.rs +++ b/crates/specialized/rtx-fea/src/analysis/mod.rs @@ -7,6 +7,7 @@ //! coordinating all lower-level components into complete workflows. pub mod dynamic_analysis; +pub mod flag3d; pub mod modal_analysis; pub mod nonlinear_analysis; pub mod nonlinear_dynamic; diff --git a/crates/specialized/rtx-fea/tests/flag3d_structure.rs b/crates/specialized/rtx-fea/tests/flag3d_structure.rs new file mode 100644 index 0000000..38ea1fd --- /dev/null +++ b/crates/specialized/rtx-fea/tests/flag3d_structure.rs @@ -0,0 +1,666 @@ +//! R8-b: the Turek–Hron flag as a 3-D Hex20 total-Lagrangian SVK solid +//! ([`rtx_fea::analysis::flag3d`]). +//! +//! Suite (fast, run by default): +//! +//! 1. `flag3d_mesh_mass_and_surface_forces` — node/element counts, the +//! consistent mass sums to `ρ V`, and the face integrator's totals +//! (uniform traction on the top face = `t · L · span`; a uniform +//! pressure on bottom + top cancels; a rigid translation of the +//! configuration changes nothing). +//! 2. `plane_strain_3d_reproduces_the_2d_csm1` — CSM1 (static, gravity) +//! with `u_z = 0` everywhere reproduces the 2-D 35×2 Quad8 plane-strain +//! model to rounding (same Newton, same banded LU). +//! 3. `plane_strain_3d_reproduces_the_2d_csm3_start` — the first 60 +//! Newmark steps of CSM3 agree with the 2-D stepper to rounding. +//! +//! Instruments (`#[ignore]`, env-driven, write under `FLAG3D_OUT`): +//! +//! * `flag3d_csm1_table` — CSM1 tip displacement per configuration. +//! * `flag3d_csm3_march` — the full CSM3 oscillation, CSV of point A and +//! of the tip's lateral corners. +//! * `flag3d_modes_dump` — the linearised operators (TL tangent at u = 0, +//! consistent mass) on the free DOFs, for an outside eigen-solve. + +use std::io::Write as _; + +use nalgebra::{DVector, Vector3}; +use rtx_fea::analysis::flag3d::{Flag3d, Flag3dSpec, FlagSide, LateralFaces}; +use rtx_fea::analysis::{ + ConvergenceCriteria, DynamicState, NonlinearDynamicAnalysis, NonlinearDynamicStepper, +}; +use rtx_fea::assembly::dof_mapping::DofComponent; +use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType}; +use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction}; +use rtx_fea::elements::total_lagrangian::{internal_force_and_tangent, saint_venant_kirchhoff}; +use rtx_fea::elements::{ElementMatrixComputer, StandardFiniteElement}; +use rtx_fea::materials::{LinearElastic, Material as _}; +use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId}; + +const E_MOD: f64 = 1.4e6; +const NU: f64 = 0.4; +/// CSM1/CSM3 density and gravity (FSI2's structure is ρ_s = 1e4). +const RHO_CSM: f64 = 1000.0; +const G: f64 = 2.0; + +fn env_str(name: &str, default: &str) -> String { + std::env::var(name).unwrap_or_else(|_| default.to_string()) +} + +fn env_num(name: &str, default: f64) -> f64 { + std::env::var(name) + .map(|v| v.parse().expect(name)) + .unwrap_or(default) +} + +/// The 2-D flag, exactly as the FSI2 harness builds it. +fn quad8_flag(nx: usize, ny: usize) -> Mesh { + let (x0, x1, y0, y1) = (0.25, 0.6, 0.19, 0.21); + let mut mesh = Mesh::new(2).unwrap(); + let (lx, ly) = (2 * nx + 1, 2 * ny + 1); + let mut grid = vec![vec![None; ly]; lx]; + for (i, column) in grid.iter_mut().enumerate() { + for (j, slot) in column.iter_mut().enumerate() { + if i % 2 == 1 && j % 2 == 1 { + continue; + } + let x = x0 + (x1 - x0) * i as f64 / (2 * nx) as f64; + let y = y0 + (y1 - y0) * j as f64 / (2 * ny) as f64; + *slot = Some(mesh.add_node(Node::new_2d(x, y))); + } + } + for i in 0..nx { + for j in 0..ny { + let (a, b) = (2 * i, 2 * j); + let nodes = vec![ + grid[a][b].unwrap(), + grid[a + 2][b].unwrap(), + grid[a + 2][b + 2].unwrap(), + grid[a][b + 2].unwrap(), + grid[a + 1][b].unwrap(), + grid[a + 2][b + 1].unwrap(), + grid[a + 1][b + 2].unwrap(), + grid[a][b + 1].unwrap(), + ]; + mesh.add_element(Element::new(ElementType::Quad8, nodes, MaterialId(0)).unwrap()) + .unwrap(); + } + } + mesh +} + +fn clamp_2d(mesh: &Mesh) -> BoundaryConditionSet { + let clamped: Vec = mesh + .nodes + .iter() + .filter(|(_, node)| (node.position().x - 0.25).abs() < 1e-12) + .map(|(&id, _)| id) + .collect(); + let mut set = BoundaryConditionSet::new(); + for component in [DofComponent::DisplacementX, DofComponent::DisplacementY] { + set.add_condition(BoundaryCondition::Dirichlet(DirichletBC { + nodes: clamped.clone(), + components: vec![component], + condition_type: DirichletType::Spatial(SpatialFunction(Box::new(|_| 0.0))), + time_range: None, + ramping_factor: 1.0, + gradual_enforcement: false, + })); + } + set +} + +fn point_2d(mesh: &Mesh, x: f64, y: f64) -> NodeId { + mesh.nodes + .iter() + .find(|(_, n)| (n.position().x - x).abs() < 1e-12 && (n.position().y - y).abs() < 1e-12) + .map(|(&id, _)| id) + .unwrap() +} + +/// Consistent gravity nodal forces `∫ N_a ρ g dV` (per unit depth in 2-D). +fn gravity_forces(mesh: &Mesh, rho: f64, g: f64) -> Vec<(NodeId, Vector3)> { + let dim = mesh.spatial_dimension; + let mut acc: std::collections::BTreeMap> = Default::default(); + for element in mesh.elements.values() { + let coords: Vec> = element + .nodes + .iter() + .map(|id| mesh.get_node(*id).unwrap().position()) + .collect(); + let fe = StandardFiniteElement::new(element.element_type, coords.clone()); + let f = ElementMatrixComputer::compute_body_force_vector( + &fe, + &coords, + &|_| Vector3::new(0.0, -rho * g, 0.0), + None, + ) + .unwrap(); + for (a, id) in element.nodes.iter().enumerate() { + let e = acc.entry(*id).or_insert_with(Vector3::zeros); + for c in 0..dim { + e[c] += f[a * dim + c]; + } + } + } + acc.into_iter().collect() +} + +/// Newton to rounding: the static comparisons are otherwise limited by +/// the default 1e-6 stopping rule (whose force scale differs between the +/// 2-D per-unit-depth and the 3-D per-span loads). +fn static_criteria() -> ConvergenceCriteria { + ConvergenceCriteria { + force_tolerance: 1e-12, + displacement_tolerance: 1e-14, + max_iterations: 60, + ..ConvergenceCriteria::default() + } +} + +/// Static equilibrium through the dynamic stepper: one "Newmark step" of +/// `dt = 1e4 s` from `u = v = a = 0` is Newton on +/// `f_int(u) + M u/(β Δt²) = F` — the static problem up to a mass term +/// 1e-9 of the stiffness. `load_steps` ramps the load, each step starting +/// from the previous equilibrium with zero velocity and acceleration. +fn static_solve<'a>( + analysis: &'a NonlinearDynamicAnalysis, + forces: &[(NodeId, Vector3)], + load_steps: usize, +) -> (DVector, NonlinearDynamicStepper<'a>, usize) { + let mut stepper = analysis.stepper().unwrap(); + let n = stepper.rest_state().unwrap().displacement.len(); + let mut u = DVector::zeros(n); + let mut iterations = 0; + for s in 1..=load_steps { + let scale = s as f64 / load_steps as f64; + let scaled: Vec<_> = forces.iter().map(|(id, f)| (*id, f * scale)).collect(); + stepper.set_nodal_forces(&scaled); + let state = DynamicState { + displacement: u.clone(), + velocity: DVector::zeros(n), + acceleration: DVector::zeros(n), + }; + let (next, it) = stepper.step(&state).unwrap(); + iterations += it; + u = next.displacement; + } + (u, stepper, iterations) +} + +fn static_2d_csm1(nx: usize, ny: usize) -> (f64, f64, usize) { + let mesh = quad8_flag(nx, ny); + let a = point_2d(&mesh, 0.6, 0.2); + let forces = gravity_forces(&mesh, RHO_CSM, G); + let analysis = NonlinearDynamicAnalysis::new( + mesh.clone(), + Flag3d::materials(E_MOD, NU, RHO_CSM), + clamp_2d(&mesh), + 1e4, + 1, + Default::default(), + ) + .with_total_lagrangian() + .with_convergence_criteria(static_criteria()); + let (u, stepper, it) = static_solve(&analysis, &forces, 5); + let d = stepper.node_dofs(a); + (u[d[0]], u[d[1]], it) +} + +struct Tip3d { + a: Vector3, + /// Point A's line at the two lateral faces (z0, z1). + side_low: Vector3, + side_high: Vector3, + iterations: usize, + dofs: usize, +} + +fn static_3d_csm1(spec: Flag3dSpec, lateral: LateralFaces, load_steps: usize) -> Tip3d { + let flag = Flag3d::build(spec).unwrap(); + let forces = gravity_forces(&flag.mesh, RHO_CSM, G); + let analysis = flag + .dynamic_analysis(E_MOD, NU, RHO_CSM, lateral, 1e4, 1, 0.5) + .with_convergence_criteria(static_criteria()); + let (u, stepper, iterations) = static_solve(&analysis, &forces, load_steps); + let read = |id: NodeId| { + let d = stepper.node_dofs(id); + Vector3::new(u[d[0]], u[d[1]], u[d[2]]) + }; + let ym = 0.5 * (spec.y0 + spec.y1); + Tip3d { + a: read(flag.point_a()), + side_low: read(flag.nearest_node(Vector3::new(spec.x1, ym, spec.z0))), + side_high: read(flag.nearest_node(Vector3::new(spec.x1, ym, spec.z1))), + iterations, + dofs: u.len(), + } +} + +fn rel(a: f64, b: f64) -> f64 { + ((a - b) / b).abs() +} + +#[test] +fn flag3d_mesh_mass_and_surface_forces() { + let spec = Flag3dSpec::turek_hron(0.1, -0.05, 7, 2, 3); + let flag = Flag3d::build(spec).unwrap(); + // Serendipity lattice: points with at most one odd index. + let [di, dj, dk] = flag.lattice_dims(); + let mut expected = 0; + for i in 0..di { + for j in 0..dj { + for k in 0..dk { + if (i % 2) + (j % 2) + (k % 2) <= 1 { + expected += 1; + } + } + } + } + assert_eq!(flag.mesh.nodes.len(), expected); + assert_eq!(flag.mesh.elements.len(), 7 * 2 * 3); + // Consistent mass sums to ρ V; every Jacobian is positive. + let mut mass = 0.0; + for element in flag.mesh.elements.values() { + let coords: Vec> = element + .nodes + .iter() + .map(|id| flag.mesh.get_node(*id).unwrap().position()) + .collect(); + let fe = StandardFiniteElement::new(element.element_type, coords.clone()); + let m = + ElementMatrixComputer::compute_consistent_mass_matrix(&fe, &coords, 1e4, None).unwrap(); + mass += m.matrix.sum(); + } + let volume = 0.35 * 0.02 * 0.1; + assert!( + rel(mass, 1e4 * volume) < 1e-12, + "mass {mass} vs {}", + 1e4 * volume + ); + + // Uniform traction on the top face: total = t · L · span. + let top = flag.surface_faces(&[FlagSide::Top]); + let t0 = Vector3::new(3.0, -2.0, 0.5); + let total: Vector3 = flag + .face_nodal_forces(&top, None, &|_, _| t0) + .iter() + .map(|(_, f)| f) + .sum(); + assert!((total - t0 * (0.35 * 0.1)).norm() < 1e-12, "{total:?}"); + // Normals point out: a pressure p on the top pushes down, on the tip + // pushes −x, on the side faces pushes inward; bottom + top cancel. + let p = 7.0; + let pressure = |_: Vector3, n: Vector3| -p * n; + let sum = |sides: &[FlagSide]| -> Vector3 { + flag.face_nodal_forces(&flag.surface_faces(sides), None, &pressure) + .iter() + .map(|(_, f)| f) + .sum() + }; + assert!((sum(&[FlagSide::Top]) - Vector3::new(0.0, -p * 0.035, 0.0)).norm() < 1e-12); + assert!((sum(&[FlagSide::Tip]) - Vector3::new(-p * 0.002, 0.0, 0.0)).norm() < 1e-12); + assert!((sum(&[FlagSide::SideHigh]) - Vector3::new(0.0, 0.0, -p * 0.007)).norm() < 1e-12); + assert!((sum(&[FlagSide::SideLow]) - Vector3::new(0.0, 0.0, p * 0.007)).norm() < 1e-12); + assert!(sum(&[FlagSide::Bottom, FlagSide::Top]).norm() < 1e-12); + // All five wetted faces + the root would close; without the root the + // pressure resultant is the root's missing +x share. + let wetted: Vector3 = flag + .face_nodal_forces(&flag.wetted_faces(), None, &pressure) + .iter() + .map(|(_, f)| f) + .sum(); + assert!((wetted - Vector3::new(-p * 0.002, 0.0, 0.0)).norm() < 1e-12); + // A rigid translation of the configuration changes nothing. + let analysis = flag.dynamic_analysis(1.4e6, 0.4, 1e4, LateralFaces::Free, 1e-3, 1, 0.5); + let stepper = analysis.stepper().unwrap(); + let mut u = DVector::zeros(3 * flag.mesh.nodes.len()); + for id in flag.mesh.nodes.keys() { + let d = stepper.node_dofs(*id); + u[d[0]] = 0.01; + u[d[1]] = -0.03; + u[d[2]] = 0.02; + } + let dofs = |id: NodeId| stepper.node_dofs(id); + let moved = flag.face_nodal_forces(&top, Some((&u, &dofs)), &|_, n| -p * n); + let still = flag.face_nodal_forces(&top, None, &|_, n| -p * n); + for ((ia, fa), (ib, fb)) in moved.iter().zip(&still) { + assert_eq!(ia, ib); + assert!((fa - fb).norm() < 1e-14); + } +} + +#[test] +fn plane_strain_3d_reproduces_the_2d_csm1() { + let (ux2, uy2, it2) = static_2d_csm1(35, 2); + let tip = static_3d_csm1( + Flag3dSpec::turek_hron(0.05, 0.0, 35, 2, 1), + LateralFaces::PlaneStrain, + 5, + ); + println!( + " CSM1 35x2: 2-D Quad8 u(A) = ({ux2:.9e}, {uy2:.9e}) [{it2} Newton]; 3-D Hex20 \ + 35x2x1 plane strain u(A) = ({:.9e}, {:.9e}, {:.2e}) [{} Newton, {} DOFs]; \ + reference (−7.18777e-3, −66.1029e-3)", + tip.a.x, tip.a.y, tip.a.z, tip.iterations, tip.dofs + ); + assert!(rel(tip.a.x, ux2) < 1e-8, "ux {} vs 2-D {ux2}", tip.a.x); + assert!(rel(tip.a.y, uy2) < 1e-8, "uy {} vs 2-D {uy2}", tip.a.y); + assert!(tip.a.z.abs() < 1e-15); + // Span-uniform: both lateral faces carry the mid-span value. + assert!((tip.side_low - tip.a).norm() < 1e-9 * tip.a.norm()); + assert!((tip.side_high - tip.a).norm() < 1e-9 * tip.a.norm()); + // And the 2-D model is the one pinned against FEATFLOW (1% short in + // u_y at 35x2, total_lagrangian_svk.rs). + assert!(rel(uy2, -66.1029e-3) < 0.02 && rel(ux2, -7.18777e-3) < 0.04); +} + +/// The free-lateral-face path, pinned: a narrow strip (span 0.02 = the +/// thickness) under CSM1 gravity. Measured with `flag3d_csm1_table` +/// (R8-b, 2026-09-25): u_y(A) = −76.340e-3 — softer than plane strain +/// (−65.141e-3) because the free faces relax the spanwise stress. +#[test] +fn free_lateral_faces_csm1_strip_pin() { + let tip = static_3d_csm1( + Flag3dSpec::turek_hron(0.02, -0.01, 35, 2, 1), + LateralFaces::Free, + 5, + ); + println!( + " CSM1 35x2x1 span 0.02 free faces: u(A) = ({:.6e}, {:.6e}, {:.2e})", + tip.a.x, tip.a.y, tip.a.z + ); + assert!(rel(tip.a.y, -76.340_06e-3) < 1e-5, "uy {}", tip.a.y); + assert!(rel(tip.a.x, -9.680_464e-3) < 1e-5, "ux {}", tip.a.x); + assert!( + tip.a.z.abs() < 1e-12, + "mid-span must not move in z: {}", + tip.a.z + ); + assert!( + (tip.side_low.y - tip.side_high.y).abs() < 1e-12, + "span symmetry" + ); +} + +fn csm3_2d(dt: f64) -> (NonlinearDynamicAnalysis, NodeId) { + let mesh = quad8_flag(35, 2); + let a = point_2d(&mesh, 0.6, 0.2); + let mut analysis = NonlinearDynamicAnalysis::new( + mesh.clone(), + Flag3d::materials(E_MOD, NU, RHO_CSM), + clamp_2d(&mesh), + dt, + 1, + Default::default(), + ) + .with_total_lagrangian(); + analysis.set_body_force(|_| Vector3::new(0.0, -RHO_CSM * G, 0.0)); + (analysis, a) +} + +#[test] +fn plane_strain_3d_reproduces_the_2d_csm3_start() { + let dt = 0.005; + let steps = 60; + let (a2d, node2) = csm3_2d(dt); + let mut s2 = a2d.stepper().unwrap(); + let flag = Flag3d::build(Flag3dSpec::turek_hron(0.05, 0.0, 35, 2, 1)).unwrap(); + let mut a3d = flag.dynamic_analysis(E_MOD, NU, RHO_CSM, LateralFaces::PlaneStrain, dt, 1, 0.5); + a3d.set_body_force(|_| Vector3::new(0.0, -RHO_CSM * G, 0.0)); + let mut s3 = a3d.stepper().unwrap(); + let node3 = flag.point_a(); + let (d2, d3) = (s2.node_dofs(node2), s3.node_dofs(node3)); + let mut st2 = s2.rest_state().unwrap(); + let mut st3 = s3.rest_state().unwrap(); + let mut worst: f64 = 0.0; + let mut peak: f64 = 0.0; + for _ in 0..steps { + st2 = s2.step(&st2).unwrap().0; + st3 = s3.step(&st3).unwrap().0; + for c in 0..2 { + worst = worst.max((st2.displacement[d2[c]] - st3.displacement[d3[c]]).abs()); + peak = peak.max(st2.displacement[d2[c]].abs()); + } + } + println!( + " CSM3 first {steps} steps (t = {:.2} s): max |u_3D − u_2D| at A {worst:.3e} m, \ + peak |u| {peak:.3e} m", + steps as f64 * dt + ); + assert!(peak > 1e-2, "the flag must have moved: {peak}"); + assert!( + worst < 1e-8 * peak, + "3-D plane strain departs from 2-D: {worst:.3e}" + ); +} + +// --------------------------------------------------------------------------- +// Instruments +// --------------------------------------------------------------------------- + +/// `NXxNYxNZ:span:free|ps` entries, comma-separated. +fn parse_configs(spec: &str) -> Vec<(usize, usize, usize, f64, LateralFaces)> { + spec.split(',') + .map(|entry| { + let mut parts = entry.trim().split(':'); + let mesh = parts.next().unwrap(); + let span: f64 = parts.next().unwrap().parse().unwrap(); + let lateral = match parts.next().unwrap() { + "free" => LateralFaces::Free, + "ps" => LateralFaces::PlaneStrain, + other => panic!("lateral {other}"), + }; + let n: Vec = mesh.split('x').map(|t| t.parse().unwrap()).collect(); + (n[0], n[1], n[2], span, lateral) + }) + .collect() +} + +fn tag(nx: usize, ny: usize, nz: usize, span: f64, lateral: LateralFaces) -> String { + let l = if lateral == LateralFaces::Free { + "free" + } else { + "ps" + }; + format!("{nx}x{ny}x{nz}_s{span}_{l}") +} + +#[test] +#[ignore = "instrument: CSM1 tip displacement per configuration"] +fn flag3d_csm1_table() { + let out = env_str("FLAG3D_OUT", "."); + let configs = parse_configs(&env_str( + "FLAG3D_CONFIGS", + "35x2x1:0.05:ps,35x2x4:0.41:free", + )); + let steps = env_num("FLAG3D_LOAD_STEPS", 5.0) as usize; + let mut table = std::fs::OpenOptions::new() + .create(true) + .append(true) + .open(format!("{out}/csm1_table.txt")) + .unwrap(); + for (nx, ny, nz, span, lateral) in configs { + let start = std::time::Instant::now(); + let tip = static_3d_csm1( + Flag3dSpec::turek_hron(span, -0.5 * span, nx, ny, nz), + lateral, + steps, + ); + let line = format!( + "CSM1 {} dofs {} newton {} | A ux {:.6e} uy {:.6e} uz {:.3e} | side_low uy {:.6e} \ + side_high uy {:.6e} | {:.1} s | ref ux -7.18777e-3 uy -66.1029e-3", + tag(nx, ny, nz, span, lateral), + tip.dofs, + tip.iterations, + tip.a.x, + tip.a.y, + tip.a.z, + tip.side_low.y, + tip.side_high.y, + start.elapsed().as_secs_f64() + ); + println!("{line}"); + writeln!(table, "{line}").unwrap(); + } + if env_str("FLAG3D_2D", "1") == "1" { + for (nx, ny) in [(35, 2), (70, 4)] { + let (ux, uy, it) = static_2d_csm1(nx, ny); + let line = + format!("CSM1 2-D {nx}x{ny} Quad8 | A ux {ux:.6e} uy {uy:.6e} [{it} Newton]"); + println!("{line}"); + writeln!(table, "{line}").unwrap(); + } + } +} + +#[test] +#[ignore = "instrument: the CSM3 oscillation on the 3-D flag"] +fn flag3d_csm3_march() { + let out = env_str("FLAG3D_OUT", "."); + let configs = parse_configs(&env_str("FLAG3D_CONFIGS", "35x2x1:0.05:ps")); + let dt = env_num("FLAG3D_DT", 0.005); + let steps = env_num("FLAG3D_STEPS", 2000.0) as usize; + for (nx, ny, nz, span, lateral) in configs { + let spec = Flag3dSpec::turek_hron(span, -0.5 * span, nx, ny, nz); + let flag = Flag3d::build(spec).unwrap(); + let mut analysis = flag.dynamic_analysis(E_MOD, NU, RHO_CSM, lateral, dt, steps, 0.5); + analysis.set_body_force(|_| Vector3::new(0.0, -RHO_CSM * G, 0.0)); + let mut stepper = analysis.stepper().unwrap(); + let ym = 0.5 * (spec.y0 + spec.y1); + let probes = [ + flag.point_a(), + flag.nearest_node(Vector3::new(spec.x1, ym, spec.z0)), + flag.nearest_node(Vector3::new(spec.x1, ym, spec.z1)), + ]; + let dofs: Vec> = probes.iter().map(|p| stepper.node_dofs(*p)).collect(); + let name = tag(nx, ny, nz, span, lateral); + let path = format!("{out}/csm3_{name}_dt{dt}.csv"); + let mut csv = std::fs::File::create(&path).unwrap(); + writeln!(csv, "t,ax,ay,az,low_y,high_y,low_z,high_z,newton").unwrap(); + let mut state = stepper.rest_state().unwrap(); + let start = std::time::Instant::now(); + let mut total = 0usize; + for k in 1..=steps { + let (next, it) = stepper.step(&state).unwrap(); + state = next; + total += it; + let u = &state.displacement; + writeln!( + csv, + "{:.6e},{:.12e},{:.12e},{:.12e},{:.12e},{:.12e},{:.12e},{:.12e},{it}", + k as f64 * dt, + u[dofs[0][0]], + u[dofs[0][1]], + u[dofs[0][2]], + u[dofs[1][1]], + u[dofs[2][1]], + u[dofs[1][2]], + u[dofs[2][2]] + ) + .unwrap(); + } + println!( + "CSM3 {name} dt {dt}: {steps} steps, {total} Newton, rescues {:?}, {:.0} s → {path}", + stepper.rescue_counts(), + start.elapsed().as_secs_f64() + ); + } +} + +#[test] +#[ignore = "instrument: the 3-D flag's K and M for an outside eigen-solve"] +fn flag3d_modes_dump() { + let out = env_str("FLAG3D_OUT", "."); + let configs = parse_configs(&env_str("FLAG3D_CONFIGS", "35x2x1:0.05:ps")); + let rho = env_num("FLAG3D_RHO", 1e4); + let (lambda, mu) = LinearElastic::new(E_MOD, NU) + .with_density(rho) + .properties() + .lame_parameters(); + let constitutive = saint_venant_kirchhoff(lambda, mu, 3); + for (nx, ny, nz, span, lateral) in configs { + let spec = Flag3dSpec::turek_hron(span, -0.5 * span, nx, ny, nz); + let flag = Flag3d::build(spec).unwrap(); + let mut ids: Vec = flag.mesh.nodes.keys().copied().collect(); + ids.sort(); + let mut free = std::collections::HashMap::new(); + let mut dof_lines = Vec::new(); + for id in &ids { + let p = flag.mesh.get_node(*id).unwrap().position(); + if (p.x - spec.x0).abs() < 1e-12 { + continue; + } + let comps = if lateral == LateralFaces::PlaneStrain { + 2 + } else { + 3 + }; + for c in 0..comps { + free.insert((*id, c), dof_lines.len()); + dof_lines.push(format!( + "{} {} {:.12e} {:.12e} {:.12e} {c}", + dof_lines.len(), + id.0, + p.x, + p.y, + p.z + )); + } + } + let mut k_trip: std::collections::BTreeMap<(usize, usize), f64> = Default::default(); + let mut m_trip: std::collections::BTreeMap<(usize, usize), f64> = Default::default(); + for element in flag.mesh.elements.values() { + let coords: Vec> = element + .nodes + .iter() + .map(|id| flag.mesh.get_node(*id).unwrap().position()) + .collect(); + let fe = StandardFiniteElement::new(element.element_type, coords.clone()); + let zero = DVector::zeros(3 * element.nodes.len()); + let (_, k_e) = + internal_force_and_tangent(&fe, &coords, &zero, constitutive.as_ref(), None) + .unwrap(); + let m_s = + ElementMatrixComputer::compute_consistent_mass_matrix(&fe, &coords, rho, None) + .unwrap(); + let local: Vec> = element + .nodes + .iter() + .flat_map(|n| (0..3).map(move |c| (*n, c))) + .map(|key| free.get(&key).copied()) + .collect(); + for (a, ga) in local.iter().enumerate() { + let Some(ga) = ga else { continue }; + for (b, gb) in local.iter().enumerate() { + let Some(gb) = gb else { continue }; + *k_trip.entry((*ga, *gb)).or_default() += k_e[(a, b)]; + if a % 3 == b % 3 { + *m_trip.entry((*ga, *gb)).or_default() += m_s.matrix[(a / 3, b / 3)]; + } + } + } + } + let name = tag(nx, ny, nz, span, lateral); + let write = |kind: &str, trip: &std::collections::BTreeMap<(usize, usize), f64>| { + let mut f = std::fs::File::create(format!("{out}/{kind}_{name}.coo")).unwrap(); + for ((i, j), v) in trip { + if *v != 0.0 { + writeln!(f, "{i} {j} {v:.17e}").unwrap(); + } + } + }; + write("k", &k_trip); + write("m", &m_trip); + std::fs::write( + format!("{out}/dofs_{name}.txt"), + dof_lines.join("\n") + "\n", + ) + .unwrap(); + println!( + "modes dump {name}: {} free DOFs, K nnz {}", + dof_lines.len(), + k_trip.len() + ); + } +}