diff --git a/ROADMAP.md b/ROADMAP.md index 0f9a3ce..cc1cab7 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -33,10 +33,10 @@ pressure) is a degree of freedom of the existing DOF framework, not a separate s | 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 — `fieldcraft-fem` kernels, `fieldcraft-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; 1D elements I: rod2/rod3 trusses, point masses, rotary inertia, springs and dashpots in both solvers and all three decks, verified by the method of joints, the three-bar truss, spring–mass, two-mass chain and torsional pendulum frequencies, and damped decay; 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); shells and beams (done: quad4 MITC4 and tri3 MITC3 flat shells with a Hughes–Brezzi drilling rotation, two-node Timoshenko beams with profile properties and orientations, six degrees of freedom at shell and beam nodes shared with solids and trusses, section forces and surface stresses, verified by shell patch tests, Navier's plate, cantilever and portal-frame beams, and the Scordelis–Lo and pinched-cylinder benchmarks); 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 (`fieldcraft-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 (trivariate TFI, curved faces) and sweep meshing (multi-level extrusions with imprinted cross-sections, revolutions, one-to-one sweeps through linking faces; interval assignment, hex8/hex20 on the geometry, B-rep face tags, `mesh.hex.generate`/`sweep`/`mapped`: done, 16 of 21 corpus parts all-hex); conforming interfaces between cells (done: `mesh.hex.blocks` meshes a part split with SolidCraft's Split Body into mapped blocks as one conforming mesh — composite faces for seams and tangent faces, interval assignment over all cells, shared faces meshed once, grading off holes and fillets, free-internal-face check; the filleted bracket in 30 cells, hex20 within 0.1 % of a fine tet10 tip deflection); sweepable cells in a conforming mesh and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | | M6 Automatic hex | Octree/grid all-hex for any part (first version: done — `fieldcraft-mesh::octree`, `mesh.hex.octree`: 3-refinement grid with conforming templates refined by curvature, thickness and narrow faces, a core cleared to a manifold boundary, feature map onto it, padding and a pillow per patch capturing every sharp curve, foldover-free fitting, certified validity, B-rep face tags; the corpus hex column and the filleted-bracket oracle case; next: transitions for concave refined regions so thin and curved parts stop over-refining); hex-dominant by tet recombination; polycube for blocky parts; layers through thin walls; mesh editing and orphan-mesh import | started | -| M7 Assemblies, contact and friction | Needs M15. Part instances from SolidCraft occurrences; tie between non-matching meshes (mortar); implicit contact: surface-to-surface mortar and node-to-surface, penalty and augmented Lagrange, self-contact, automatic pair detection with a bounding-volume search shared with explicit contact; friction: Coulomb with stick–slip, static/kinetic, anisotropic, pressure-, rate- and temperature-dependent coefficients; contact output (pressure, slip, status); thermal contact conductance and frictional heating; wear (Archard); bolt pretension, connectors with behaviour (elastic, damped, stops, locks); amplitudes, multiple steps, restart; nonlinear steps in all three decks. Verified by Hertz (sphere and cylinder), the Cattaneo–Mindlin partial-slip problem, the sliding block, NAFEMS contact benchmarks | planned | +| M7 Assemblies, contact and friction | Needs M15. Part instances from SolidCraft occurrences; tie between non-matching meshes (mortar); implicit contact: surface-to-surface mortar and node-to-surface, penalty and augmented Lagrange, self-contact, automatic pair detection (the model's contact definitions, `interaction.contact.*` with `find`, the deck mapping, `fieldcraft_fem::contact`'s closest-point projection and bucket grid, and the explicit solver's segmentation and friction return map are in place to reuse; the implicit solver refuses contact by name until then); friction: Coulomb with stick–slip, static/kinetic, anisotropic, pressure-, rate- and temperature-dependent coefficients; contact output (pressure, slip, status); thermal contact conductance and frictional heating; wear (Archard); bolt pretension, connectors with behaviour (elastic, damped, stops, locks); amplitudes, multiple steps, restart; nonlinear steps in all three decks. Verified by Hertz (sphere and cylinder), the Cattaneo–Mindlin partial-slip problem, the sliding block, NAFEMS contact benchmarks | planned | | M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; hybrid (mixed u–p) hex8/hex20/tet10 for nearly incompressible materials (done); heat transfer (solids done: steady and θ-method transient conduction with flux, convection and sources, thermal strain in static steps (solids and 2D) and explicit dynamics (with plasticity), coupled temperature–displacement, `solve.heat`; deck output, 2D/shell/beam conduction planned) | started | | M9 Post-processing and deck readers | Results readers for public formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`); XY and path plots, cuts, display groups, free-body forces; deck readers that keep unknown keywords; mesh convergence studies | planned | -| M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `fieldcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (X2: done — updated Lagrangian elements, Jaumann/Hughes–Winget and Green–Naghdi/Flanagan–Taylor rates, J2 bilinear/table plasticity, Johnson–Cook strength with heating and damage, erosion; materials gain `plastic`, `damage`, `erosion` with `model.material.plasticity`/`damage`; Dienes shear, rigid rotation, J2 uniaxial and Taylor-bar tests); contact (global search, penalty with friction, rigid walls); runs on hex and tet meshes | started | +| M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `fieldcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (X2: done — updated Lagrangian elements, Jaumann/Hughes–Winget and Green–Naghdi/Flanagan–Taylor rates, J2 bilinear/table plasticity, Johnson–Cook strength with heating and damage, erosion; materials gain `plastic`, `damage`, `erosion` with `model.material.plasticity`/`damage`; Dienes shear, rigid rotation, J2 uniaxial and Taylor-bar tests); contact (X3: done — `interaction.contact.pair`/`general`/`find`, segment-to-segment penalty contact both ways with polygon clipping (passes the contact patch test on non-matching meshes), pairs, general (all exterior faces) and self contact, bucket-grid search, deterministic on any thread count, a penalty from the masses and the stable increment that never shortens it, Coulomb friction with exponential decay by a return map, friction on rigid walls, contact energy and friction dissipation in the balance, initial penetrations adjusted or refused; Abaqus `*CONTACT PAIR`/`*CONTACT`/`*FRICTION` and LS-DYNA automatic contacts in the decks; verified by colliding bars, sliding blocks, the patch test and Hertz's sphere impact); hyperelastic materials on the shared material point and mass scaling next; runs on hex and tet meshes | started | | M11 Equations of state and shells | Polynomial, ideal-gas, Mie–Grüneisen, JWL with programmed burn, tabulated; energy coupling; shock verification (Rankine–Hugoniot, Sod, Noh, Sedov); Belytschko–Lin–Tsay shells with self-contact on quad meshes | planned | | M12 CFD: lattice Boltzmann | D3Q19/D3Q27 on grids voxelised from the face-tagged surface; TRT, regularised and cumulant collision; interpolated walls, inlets and outlets; LES (Smagorinsky, WALE) with wall models; thermal LBM and conjugate heat transfer with the structure's thermal solver; non-Newtonian viscosity (power-law, Carreau–Yasuda, Cross, Herschel–Bulkley, Bingham) and temperature-dependent viscosity; porous media (Brinkman, Darcy–Forchheimer); multiphase and free surface (M13); validation by the lid-driven cavity, backward-facing step, cylinder wake (Strouhal) and turbulent channel; octree refinement shared with the hex mesher; wgpu backend (deterministic per device, checked against the CPU) (CPU core: done — `fieldcraft-lbm` with D3Q19 TRT and Guo forcing, halfway bounce-back and moving walls, slip and periodic faces, Zou–He velocity and pressure faces, exact voxelisation of obstacle bodies, momentum-exchange forces per face; fluid materials, fluid domains and fluid steps in the model; `model.fluid_domain.*`, `solve.fluid`; verified by Poiseuille channel and pipe, Couette, Taylor–Green decay, sphere-array Stokes drag and mass conservation; Schäfer–Turek 2D-1 as a benchmark) | started | | M13 Fluid–structure interaction | Immersed-boundary coupling of the lattice Boltzmann fluid with the explicit or implicit structure; interface force transfer; sub-cycling; Aitken then interface quasi-Newton coupling; free-surface flow | planned | @@ -77,7 +77,7 @@ What exists and what is missing, by area; each gap is in a milestone above. | Area | Today | Missing | |---|---|---| | Solution procedures | Linear static, nonlinear static (solids with geometric and material nonlinearity: hyperelasticity, J2 plasticity), frequency (Lanczos), steady and transient heat, coupled temperature–displacement, explicit dynamics, LBM flow | Nonlinear shells and beams, arc-length, buckling (M15); implicit, harmonic, random and spectrum dynamics (M16) | -| Contact, friction, adhesion | Explicit planar rigid walls (frictionless); couplings, springs, dashpots | Implicit and explicit deformable contact, tie, friction laws, thermal contact, wear (M7); cohesive adhesion and debonding (M20) | +| Contact, friction, adhesion | Explicit deformable contact (pairs, general, self; segment to segment, penalty) with Coulomb friction and exponential decay; rigid walls with friction; couplings, springs, dashpots | Implicit contact, tie, anisotropic and pressure- or temperature-dependent friction, thermal contact, wear, eroding contact (M7, M10); cohesive adhesion and debonding (M20) | | Materials | Isotropic and orthotropic elasticity, density, thermal, J2 plasticity (bilinear, table, Johnson–Cook), Johnson–Cook damage and erosion, compressible neo-Hookean, Mooney–Rivlin and Yeoh hyperelasticity, Newtonian fluids | Elastomers beyond Yeoh (Ogden, Arruda–Boyce and the rest), viscoelasticity, viscoplasticity and creep, anisotropic and pressure-dependent plasticity, non-Newtonian viscosity, a cited material library and calibration (M17) | | Composites and anisotropy | Orthotropic constants | Orientations, layups, ply results, failure criteria, homogenisation (M18) | | Coupled fields | Thermal–structural (one way and coupled) | Piezoelectric, poroelastic, electro-thermal, acoustics (M19) | diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs index b1ebe58..841fe96 100644 --- a/crates/engine/src/cmd/define.rs +++ b/crates/engine/src/cmd/define.rs @@ -266,7 +266,7 @@ fn fluid_boundary(s: &mut Session, p: &Value) -> Result { pub const STEP_CREATE: CommandSpec = CommandSpec { id: "model.step.create", label: "New Step", - params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "nonlinear_static", "nonlinear_geometry": true, "initial_increment": 0.1, "min_increment": 1e-5, "max_increment": 1.0 (fractions of the step), "max_increments": 100, "max_iterations": 16 (per increment), "tolerances": {"force": 1e-6, "displacement": 1e-5, "energy": 1e-10}, "line_search": false, "load_path": [[0, 0], [0.5, 1], [1, 0]] (optional: the load factor against the step fraction, piecewise linear from [0, 0] to fraction 1; increments end at its points)} (loads and prescribed displacements are the load factor times their values; it ramps from 0 to 1 without a load path) | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500} | {"type": "fluid", "domain": "Flow", "end_time": 1.0 | "steps": 1000, "time_step": (default: chosen for Mach 0.087 and τ ≤ 1), "frames": 10, "magic": 0.1875 (the two-relaxation-time Λ)} | {"type": "explicit_dynamic", "end_time": 0.001, "dt_scale": 0.9, "max_increments": 10000000, "output_interval", "history_interval", "history_nodes": "set", "initial_velocities": [{"set", "velocity": [vx, vy, vz]}], "rigid_walls": [{"set", "point", "normal"}], "hourglass": {"type": "viscous|stiffness", "coefficient": 0.1}, "bulk_viscosity": {"linear": 0.06, "quadratic": 1.5}}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, + params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "nonlinear_static", "nonlinear_geometry": true, "initial_increment": 0.1, "min_increment": 1e-5, "max_increment": 1.0 (fractions of the step), "max_increments": 100, "max_iterations": 16 (per increment), "tolerances": {"force": 1e-6, "displacement": 1e-5, "energy": 1e-10}, "line_search": false, "load_path": [[0, 0], [0.5, 1], [1, 0]] (optional: the load factor against the step fraction, piecewise linear from [0, 0] to fraction 1; increments end at its points)} (loads and prescribed displacements are the load factor times their values; it ramps from 0 to 1 without a load path) | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500} | {"type": "fluid", "domain": "Flow", "end_time": 1.0 | "steps": 1000, "time_step": (default: chosen for Mach 0.087 and τ ≤ 1), "frames": 10, "magic": 0.1875 (the two-relaxation-time Λ)} | {"type": "explicit_dynamic", "end_time": 0.001, "dt_scale": 0.9, "max_increments": 10000000, "output_interval", "history_interval", "history_nodes": "set", "initial_velocities": [{"set", "velocity": [vx, vy, vz]}], "rigid_walls": [{"set", "point", "normal", "friction": {"static": 0.3, "kinetic": 0.2, "decay": 10} (Coulomb; default frictionless)}], "hourglass": {"type": "viscous|stiffness", "coefficient": 0.1}, "bulk_viscosity": {"linear": 0.06, "quadratic": 1.5}}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, kind: CommandKind::Model, run: step_create, }; diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index c4b8446..9eea63d 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -58,7 +58,7 @@ pub const FREQUENCY: CommandSpec = CommandSpec { pub const EXPLICIT: CommandSpec = CommandSpec { id: "solve.explicit", label: "Solve Explicit Dynamic Step", - params: r#"{"step": "Impact", "max_dofs": 6000000, "max_increments": 1000000000, "max_output_values": 50000000, "threads": 8} — transient response of an explicit dynamic step (the only one when `step` is left out) by central differences; keeps its histories and field frames in the session; the caps may only be lowered (the step's own increment cap applies too); `threads` as for solve.static"#, + params: r#"{"step": "Impact", "max_dofs": 6000000, "max_increments": 1000000000, "max_output_values": 50000000, "max_contact_faces": 10000000, "threads": 8} — transient response of an explicit dynamic step (the only one when `step` is left out) by central differences, with the model's contacts (interaction.contact.*) and the step's rigid walls; keeps its histories and field frames in the session and reports each contact's peak force, first and last contact, deepest penetration, friction dissipated and the nodes adjusted at the start; the caps may only be lowered (the step's own increment cap applies too); `threads` as for solve.static"#, kind: CommandKind::Model, run: run_explicit, }; @@ -66,7 +66,7 @@ pub const EXPLICIT: CommandSpec = CommandSpec { pub const HISTORY: CommandSpec = CommandSpec { id: "solve.history", label: "Result Histories", - params: r#"{"step": "Impact", "every": 1} — the time histories of a solved explicit step: time, increment, energies, momentum, constraint force and the history nodes' displacement and velocity; `every` keeps one sample in so many (by default enough to return at most 10000)"#, + params: r#"{"step": "Impact", "every": 1} — the time histories of a solved explicit step: time, increment, energies, momentum, constraint force, the normal force across each contact (in `contacts` order) and the history nodes' displacement and velocity; `every` keeps one sample in so many (by default enough to return at most 10000)"#, kind: CommandKind::Model, run: history, }; @@ -289,6 +289,7 @@ fn run_explicit(s: &mut Session, p: &Value) -> Result { max_dofs: cap(p, "max_dofs", ExplicitLimits::MAX_DOFS)?, max_increments: cap(p, "max_increments", ExplicitLimits::MAX_INCREMENTS)?, max_output_values: cap(p, "max_output_values", ExplicitLimits::MAX_OUTPUT_VALUES)?, + max_contact_faces: cap(p, "max_contact_faces", ExplicitLimits::MAX_CONTACT_FACES)?, }; let step = opt::(p, "step")?; let model = &s.model; @@ -314,6 +315,7 @@ fn explicit_summary(r: &ExplicitResult) -> Value { "frames": r.frames.len(), "history_samples": r.history.time.len(), "eroded": r.eroded, + "contacts": r.contacts, "warnings": r.warnings, }) } @@ -355,6 +357,8 @@ fn history(s: &mut Session, p: &Value) -> Result { "energies": pick(&|i| json!(h.energies.get(i))), "momentum": pick(&|i| json!(h.momentum.get(i))), "constraint_force": pick(&|i| json!(h.constraint_force.get(i))), + "contact_force": pick(&|i| json!(h.contact_force.get(i))), + "contacts": r.contacts.iter().map(|c| c.name.as_str()).collect::>(), "nodes": nodes, })) } diff --git a/crates/explicit/src/contact.rs b/crates/explicit/src/contact.rs new file mode 100644 index 0000000..924f16a --- /dev/null +++ b/crates/explicit/src/contact.rs @@ -0,0 +1,822 @@ +//! Penalty contact between element faces, with Coulomb friction. +//! +//! **Kinematics: segment to segment, both ways.** The normal contact pressure is integrated over +//! the overlap of each pair of facing faces: a secondary face and a main face it faces are +//! projected onto the secondary face's plane, each split into triangles, clipped against each +//! other (Sutherland–Hodgman), and the overlap is integrated with a three-point rule per +//! sub-triangle; at each point the gap g = (x_s − x_m) · n_m between the two faces (n_m the main +//! face's outward normal) gives a penalty pressure p = ε (−g) where they overlap, which acts on +//! the secondary face's nodes with weights N_i^s and on the main face's with −N_j^m. This is the +//! segmentation and integration of mortar contact (M. A. Puso, T. A. Laursen, *A mortar +//! segment-to-segment contact method for large deformation solid mechanics*, CMAME +//! 193:601–629, 2004) with the penalty applied at the integration points (K. A. Fischer, P. +//! Wriggers, *Frictionless 2D contact formulations for finite deformations based on the mortar +//! method*, Comput. Mech. 36:226–244, 2005): for a uniform pressure across flat faces both sides +//! receive exactly the consistent nodal forces ∫ p N dA, whatever the meshes, so the contact +//! patch test is passed — node-to-segment contact fails it on non-matching meshes (R. L. Taylor, +//! P. Papadopoulos, *On a patch test for contact problems in two dimensions*, 1991). Every +//! contact runs both ways — a pair with each surface as the secondary in turn, single-surface +//! (general and self) contact with every face against every other one that shares no node with +//! it — at half the penalty each, so neither side is privileged and no corner of either surface +//! sinks into the other. Two faces are in contact where they face each other (their normals more +//! than [`FACING`] from parallel, which keeps the faces of one body that meet at an edge apart) +//! and the gap is negative but no deeper than half the main face's element is thick (deeper, the +//! surfaces have passed through each other). Faces of eroded elements leave the contact. +//! +//! **Search.** Every increment: each surface's faces go into a bucket grid with their boxes grown +//! by their depth limit ([`fieldcraft_fem::contact::Grid`]: D. J. Benson, J. O. Hallquist, *A +//! single surface contact algorithm for the post-buckling analysis of shell structures*, CMAME +//! 78:141–163, 1990); each secondary face takes the main faces listed in the cells its box +//! overlaps. The faces are integrated in parallel, each on its own, and their forces added face by +//! face in a fixed order, so the result does not depend on how the faces are split among threads. +//! +//! **Penalty.** The penalty per area comes from the masses and the stable increment Δt: +//! ε = s γ min_i (m_i / A_i) / Δt², with m_i the lumped mass of node i and A_i the area that +//! bounds its contact stiffness, over the nodes of the contact. The stiffness entries of node i +//! in direction d sum to at most 2 ε Σ_f a_f |n_f,d| Σ_d′ |n_f,d′| over its faces f (a_f its +//! share of the face, a quarter of a quadrilateral or a third of a triangle, n_f the unit normal; +//! its weights over a face sum to one on both sides), so A_i = κ_i max_d Σ_f a_f |n_f,d| with +//! κ_i = max_f Σ_d |n_f,d| ∈ [1, √3]: the node's area on a flat surface in any orientation, and +//! no more than the largest of its areas facing one way at an edge or a corner. By Gershgorin's +//! theorem the contact then adds at most 2 s γ / Δt² to the square of any frequency (Weyl's +//! inequality adds it to the elements'); with γ = [`GAMMA`] and the +//! default increment scale (0.9 of the critical increment, ω Δt ≤ 1.8) the step stays under the +//! stability limit ω Δt = 2 (1.8² + 2γ ≤ 4) and contact never shortens the increment. The penalty +//! scale s ≤ 1 only softens it. The pressure stiffness this gives is about a fifth of the +//! smallest element's M/h (M the constrained modulus): stiffer than the element-based penalty +//! f K A²/V with its usual f = 0.1 (J. O. Hallquist, G. L. Goudreau, D. J. Benson, *Sliding +//! interfaces with contact-impact in large-scale Lagrangian computations*, CMAME 51:107–137, +//! 1985), and a uniform ε keeps a uniform pressure uniform. +//! +//! **Friction** at the secondary nodes by a regularised return map (P. Wriggers, *Computational +//! Contact Mechanics*, 2nd ed., Springer 2006, §5.2 and §10.3; J. C. Simo, T. A. Laursen, *An +//! augmented Lagrangian treatment of contact problems involving friction*, Comput. Struct. +//! 42(1):97–116, 1992): a node's normal force f_n is its share of the pressure, its tangential +//! force an elastic predictor, the previous force moved onto the current tangent plane less +//! k_t times the tangential motion of the node relative to the closest point of the main face that +//! pressed it most over the increment, t* = P t_n − k_t Δu_T with k_t the node's normal stiffness ε A_i / 2; if +//! |t*| ≤ μ f_n the node sticks (t = t*), otherwise it slips and the force returns to the cone, +//! t = μ f_n t*/|t*|, dissipating μ f_n (|t*| − μ f_n)/k_t. The reaction goes to the main face's +//! nodes with −N_j. The coefficient μ = μ_k + (μ_s − μ_k) e^{−d γ̇} follows the slip rate +//! γ̇ = |Δu_T|/Δt ([`fieldcraft_model::Friction`]). A node that leaves contact forgets its +//! tangential force. +//! +//! **Energy.** The contact forces enter the equations of motion like external forces; their work +//! (the penalty energy stored, plus what friction dissipates) is integrated with the trapezoidal +//! rule like the internal forces' (see [`crate::integrate`]), and the friction dissipation is +//! reported too. +//! +//! **Initial penetration.** Before the analysis each node of a secondary surface is projected +//! onto the faces of the main surface ([`fieldcraft_fem::contact::closest_point`]); a node behind +//! a face it lies over (within a fiftieth of the face's size), within the depth limit and with its +//! own faces facing it, is moved onto it (strain-free: the reference coordinates change before the +//! elements are built) or the analysis is refused with the node, the face and the depth, as the +//! contact's `initial_penetration` says. +//! +//! **Erosion.** The faces of eroded elements leave the contact; the faces their removal exposes +//! do not join it (eroding contact is later work). + +use fieldcraft_fem::contact::{Aabb, Grid, Projection, bounds, closest_point}; +use fieldcraft_fem::frame::{cross, dot, norm, sub}; +use fieldcraft_model::{Contact, ContactSurfaces, ElementFace, ElementKind, Friction, InitialPenetration, Model, Resolved}; +use serde::Serialize; + +use crate::elements::{State, add_at}; +use crate::system::NODE_DOFS; +use crate::{Error, Result}; + +/// The penalty's share of the stability budget: ε = s γ min(m/A) / Δt² (see the module docs). +pub const GAMMA: f64 = 0.3; + +/// A node lies over a face when its lateral offset is at most this fraction of the face's size +/// (the initial-penetration check). +pub const EDGE_TOLERANCE: f64 = 0.02; + +/// A point deeper than this fraction of the thickness of a face's element has passed through it. +pub const DEPTH_FRACTION: f64 = 0.5; + +/// Facing faces: the dot product of their unit normals is below minus this. +pub const FACING: f64 = 0.1; + +/// Faces per parallel task. +const BLOCK: usize = 64; + +/// The three-point rule on a triangle, exact for quadratics: barycentric points (of the second +/// and third corners), weights summing to one. +const TRIANGLE_RULE: [([f64; 2], f64); 3] = + [([1.0 / 6.0, 1.0 / 6.0], 1.0 / 3.0), ([2.0 / 3.0, 1.0 / 6.0], 1.0 / 3.0), ([1.0 / 6.0, 2.0 / 3.0], 1.0 / 3.0)]; + +type Vec3 = [f64; 3]; + +/// A face of a contact surface. +#[derive(Debug, Clone)] +pub struct Face { + /// The corner nodes (system indices), in the face's order; the last unused on a triangle. + pub nodes: [usize; 4], + pub count: usize, + /// The element's label, its position in [`crate::system::System::elements`] and the face + /// number. + pub element: u32, + pub index: usize, + pub face: u8, + /// The nodes of the face's element (a node of it is never in contact with the face). + pub owner: Vec, + /// The penetration beyond which a point has passed through: [`DEPTH_FRACTION`] of the + /// element's thickness under the face. + pub depth: f64, +} + +impl Face { + fn corners(&self) -> &[usize] { + self.nodes.get(..self.count).unwrap_or_default() + } +} + +/// The faces and nodes of one surface. +#[derive(Debug, Clone, Default)] +pub struct Surface { + pub faces: Vec, + /// Its nodes, sorted (system indices). + pub nodes: Vec, + /// Per node (in `nodes` order), the faces that hold it. + node_faces: Vec>, + /// Per node, the area of the surface that bounds its contact stiffness (undeformed; see + /// [`Interface::forces`]). + area: Vec, +} + +/// A node's friction state in one pass. +#[derive(Debug, Clone, Copy, Default)] +struct Slip { + /// The tangential force on the node. + force: Vec3, +} + +/// One pass: the faces of a surface (secondary) against the faces of a surface (main). +#[derive(Debug, Clone)] +struct Pass { + secondary: usize, + main: usize, + /// Per node of the secondary surface. + slip: Vec, +} + +/// A contact as the solver runs it. +#[derive(Debug, Clone)] +pub struct Interface { + pub name: String, + pub friction: Friction, + pub scale: f64, + pub surfaces: Vec, + passes: Vec, + /// Nodes moved onto the other surface before the analysis. + pub adjusted: Vec, + /// min over the contact's nodes of mass / area, once the masses are known. + mass_per_area: Option, +} + +/// A node moved onto a face before the analysis. +#[derive(Debug, Clone, Copy, PartialEq, Serialize)] +pub struct Adjustment { + pub node: u32, + pub element: u32, + pub face: u8, + /// How far it started inside. + pub depth: f64, +} + +/// What one evaluation of the contact forces found, per contact. +#[derive(Debug, Clone, Copy, Default, PartialEq)] +pub struct Evaluation { + /// The normal force across the contact (the integrated pressure, both ways). + pub normal_force: f64, + /// The deepest penetration. + pub penetration: f64, + /// Integration points in contact. + pub active: usize, + /// Friction dissipated over the increment that led here. + pub dissipation: f64, +} + +fn invalid(m: impl Into) -> Error { + Error::Invalid(m.into()) +} + +/// The volume of a solid element of corners `x` (hex8 or tet4), for the thickness under a face. +fn volume(kind: ElementKind, x: &[Vec3]) -> Option { + match kind { + ElementKind::Hex8 => <[Vec3; 8]>::try_from(x).ok().and_then(|a| fieldcraft_fem::uniform::Hex8::hex8(&a).ok()).map(|k| k.volume), + ElementKind::Tet4 => <[Vec3; 4]>::try_from(x).ok().and_then(|a| fieldcraft_fem::uniform::Tet4::tet4(&a).ok()).map(|k| k.volume), + _ => None, + } +} + +/// A face's area vector (its outward normal times its area, flat or bilinear). +fn area_vector(x: &[Vec3]) -> Vec3 { + match x { + [a, b, c] => cross(sub(*b, *a), sub(*c, *a)).map(|v| 0.5 * v), + [a, b, c, d] => cross(sub(*c, *a), sub(*d, *b)).map(|v| 0.5 * v), + _ => [0.0; 3], + } +} + +/// The largest side of a face's bounding box. +fn size(x: &[Vec3]) -> f64 { + let b = bounds(x, 0.0); + (0..3).map(|k| b[1][k] - b[0][k]).fold(0.0, f64::max) +} + +/// Where a node is: its rest position plus its displacement. +fn position(xyz: &[Vec3], u: &[f64], node: usize) -> Vec3 { + let x = xyz.get(node).copied().unwrap_or_default(); + [0, 1, 2].map(|k| x[k] + u.get(NODE_DOFS * node + k).copied().unwrap_or(0.0)) +} + +/// A node's displacement increment. +fn shift(du: &[f64], node: usize) -> Vec3 { + [0, 1, 2].map(|k| du.get(NODE_DOFS * node + k).copied().unwrap_or(0.0)) +} + +impl Surface { + /// The surface of element faces `faces`, its nodes, their faces and areas. `node` and + /// `element` give a label's index among the analysed nodes and elements. + fn new( + model: &Model, + faces: &[ElementFace], + xyz: &[Vec3], + node: &impl Fn(u32) -> Option, + element: &impl Fn(u32) -> Option, + ) -> Result { + let mut out = Surface::default(); + for f in faces { + let at = |what: &str| invalid(format!("element {} face {}: {what}", f.element, f.face)); + let el = model.element(f.element).ok_or_else(|| at("the element is missing"))?; + if !matches!(el.kind, ElementKind::Hex8 | ElementKind::Tet4) { + return Err(Error::NotSupported(format!( + "contact on the faces of {} elements (element {}); explicit contact acts on hex8 and tet4 faces", + el.kind.name(), + f.element + ))); + } + let labels = model.face_nodes(*f).ok_or_else(|| at("no such face"))?; + let corners: Vec = labels.iter().map(|l| node(*l).ok_or_else(|| at("a node is not analysed"))).collect::>()?; + let owner: Vec = el.nodes.iter().filter_map(|l| node(*l)).collect(); + let index = element(f.element).ok_or_else(|| at("the element is not analysed"))?; + let x: Vec = corners.iter().map(|n| xyz.get(*n).copied().unwrap_or_default()).collect(); + let ex: Vec = owner.iter().map(|n| xyz.get(*n).copied().unwrap_or_default()).collect(); + let area = norm(area_vector(&x)); + // Thickness V/A; a face of no area or an element of no volume falls back on its size. + let thickness = volume(el.kind, &ex).filter(|v| *v > 0.0 && area > 0.0).map_or(size(&x), |v| v / area); + let mut nodes = [usize::MAX; 4]; + for (slot, n) in nodes.iter_mut().zip(&corners) { + *slot = *n; + } + out.faces.push(Face { nodes, count: corners.len(), element: f.element, index, face: f.face, owner, depth: DEPTH_FRACTION * thickness }); + } + out.nodes = out.faces.iter().flat_map(|f| f.corners().iter().copied()).collect(); + out.nodes.sort_unstable(); + out.nodes.dedup(); + out.node_faces = vec![Vec::new(); out.nodes.len()]; + // Per node and direction, Σ_f A_f |n_f,d| / count, and the largest Σ_d |n_d| of its faces. + let mut directional = vec![[0.0; 3]; out.nodes.len()]; + let mut spread = vec![1.0f64; out.nodes.len()]; + for (i, f) in out.faces.iter().enumerate() { + let x: Vec = f.corners().iter().map(|n| xyz.get(*n).copied().unwrap_or_default()).collect(); + let a = area_vector(&x); + let l = norm(a); + let share = a.map(|c| c.abs() / f.count as f64); + let kappa = if l > 0.0 { a.iter().map(|c| c.abs()).sum::() / l } else { 1.0 }; + for n in f.corners() { + if let Ok(k) = out.nodes.binary_search(n) { + if let Some(list) = out.node_faces.get_mut(k) { + list.push(i); + } + if let (Some(d), Some(s)) = (directional.get_mut(k), spread.get_mut(k)) { + for c in 0..3 { + d[c] += share[c]; + } + *s = s.max(kappa); + } + } + } + } + out.area = directional.iter().zip(&spread).map(|(d, s)| s * d.iter().copied().fold(0.0, f64::max)).collect(); + Ok(out) + } + + /// The faces' corner positions at displacements `u`. + fn positions(&self, xyz: &[Vec3], u: &[f64]) -> Vec> { + self.faces.iter().map(|f| f.corners().iter().map(|n| position(xyz, u, *n)).collect()).collect() + } + + /// The nodes' outward normals (the sum of their faces' area vectors; not normalised). + fn node_normals(&self, faces: &[Vec]) -> Vec { + self.node_faces + .iter() + .map(|list| { + let mut n = [0.0; 3]; + for a in list.iter().filter_map(|i| faces.get(*i)).map(|x| area_vector(x)) { + for k in 0..3 { + n[k] += a[k]; + } + } + n + }) + .collect() + } +} + +/// The faces' grid for the global search at positions `positions`. +fn grid(surface: &Surface, positions: &[Vec]) -> Grid { + let boxes: Vec = surface.faces.iter().zip(positions).map(|(f, x)| bounds(x, f.depth)).collect(); + Grid::new(&boxes, 4 * boxes.len() + 64) +} + +/// The face a node at `x` (index `node`, outward normal `normal`) starts inside, among the grid's +/// candidates: over it, behind it within its depth limit, facing it; the nearest. +fn inside(surface: &Surface, positions: &[Vec], grid: &Grid, node: usize, x: Vec3, normal: Vec3) -> Option<(usize, Projection)> { + let mut best: Option<(usize, Projection)> = None; + for &c in grid.candidates(x) { + let c = c as usize; + let (Some(face), Some(fx)) = (surface.faces.get(c), positions.get(c)) else { continue }; + if face.owner.contains(&node) { + continue; + } + let Some(q) = closest_point(fx, x) else { continue }; + let over = q.lateral <= EDGE_TOLERANCE * size(fx); + let n = norm(normal); + let facing = n == 0.0 || dot(normal, q.normal) < -FACING * n; + if q.gap < 0.0 && -q.gap <= face.depth && over && facing && best.is_none_or(|b| q.distance < b.1.distance) { + best = Some((c, q)); + } + } + best +} + +/// A point of a plane in its own 2D coordinates. +type Point = [f64; 2]; + +fn cross2(o: Point, a: Point, b: Point) -> f64 { + (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0]) +} + +/// The signed area of a polygon (counter-clockwise positive). +fn signed_area(p: &[Point]) -> f64 { + let mut s = 0.0; + for (i, a) in p.iter().enumerate() { + let b = p.get((i + 1) % p.len()).copied().unwrap_or(*a); + s += a[0] * b[1] - a[1] * b[0]; + } + 0.5 * s +} + +/// `subject` clipped by the convex counter-clockwise polygon `clip` (Sutherland–Hodgman). +fn clip(subject: &[Point], clip: &[Point]) -> Vec { + let mut out: Vec = subject.to_vec(); + for (i, a) in clip.iter().enumerate() { + let Some(b) = clip.get((i + 1) % clip.len()).copied() else { break }; + let input = std::mem::take(&mut out); + let side = |p: Point| cross2(*a, b, p); + for (j, p) in input.iter().enumerate() { + let Some(q) = input.get((j + 1) % input.len()).copied() else { break }; + let (sp, sq) = (side(*p), side(q)); + if sp >= 0.0 { + out.push(*p); + } + if (sp >= 0.0) != (sq >= 0.0) { + let t = sp / (sp - sq); + out.push([p[0] + t * (q[0] - p[0]), p[1] + t * (q[1] - p[1])]); + } + } + if out.is_empty() { + break; + } + } + out +} + +/// The corners `t` of `points`, counter-clockwise. +fn counter_clockwise(t: &[usize; 3], points: &[Point]) -> Vec { + let mut out: Vec = t.iter().filter_map(|i| points.get(*i).copied()).collect(); + if signed_area(&out) < 0.0 { + out.reverse(); + } + out +} + +/// The triangles a face splits into (corner positions in its order). +fn triangles(count: usize) -> &'static [[usize; 3]] { + if count == 3 { &[[0, 1, 2]] } else { &[[0, 1, 2], [0, 2, 3]] } +} + +/// What the integration needs of each main face, worked out once per increment. +#[derive(Debug, Clone, Copy)] +struct FaceInfo { + /// Area vector and its length. + area: Vec3, + length: f64, + /// The lateral offset beyond which a point is past the face's edge. + edge: f64, + alive: bool, +} + +fn face_info(surface: &Surface, positions: &[Vec], alive: &dyn Fn(&Face) -> bool) -> Vec { + surface + .faces + .iter() + .zip(positions) + .map(|(f, x)| { + let area = area_vector(x); + FaceInfo { area, length: norm(area), edge: EDGE_TOLERANCE * size(x), alive: alive(f) } + }) + .collect() +} + +/// One secondary face against the main surface: its contributions to the nodal forces and to +/// each secondary node's normal force. +#[derive(Debug, Default)] +struct FaceForces { + /// (node, force) contributions, in order. + forces: Vec<(usize, Vec3)>, + /// (secondary node, normal force share, main face) of each integration point. + normal: Vec<(usize, f64, usize)>, + penetration: f64, + points: usize, + total: f64, +} + +/// The penalty pressure integrated over the overlaps of secondary face `f` (corners `fx`) with +/// the main faces facing it, at penalty `eps` per area. +#[allow(clippy::too_many_arguments)] +fn integrate_face( + f: &Face, + fx: &[Vec3], + alive: bool, + main: &Surface, + mx: &[Vec], + info: &[FaceInfo], + grid: &Grid, + eps: f64, + found: &mut Vec, +) -> FaceForces { + let mut out = FaceForces::default(); + if !alive { + return out; + } + let area = area_vector(fx); + let a = norm(area); + if !(a > 0.0 && a.is_finite()) { + return out; + } + let normal = area.map(|c| c / a); + let n = fx.len() as f64; + let center: Vec3 = [0, 1, 2].map(|k| fx.iter().map(|p| p[k]).sum::() / n); + let Some(first) = fx.first().copied() else { return out }; + // A basis of the face's plane, (e1, e2, normal) right-handed. + let r = sub(first, center); + let r = sub(r, normal.map(|c| c * dot(r, normal))); + let rl = norm(r); + if !(rl > 0.0) { + return out; + } + let e1 = r.map(|c| c / rl); + let e2 = cross(normal, e1); + let flat = |x: Vec3| -> Point { + let d = sub(x, center); + [dot(d, e1), dot(d, e2)] + }; + let lift = |p: Point| -> Vec3 { [0, 1, 2].map(|k| center[k] + p[0] * e1[k] + p[1] * e2[k]) }; + let own: Vec = fx.iter().map(|x| flat(*x)).collect(); + let own_triangles: Vec> = triangles(f.count).iter().map(|t| counter_clockwise(t, &own)).collect(); + grid.candidates_in(&bounds(fx, 0.0), found); + for &c in found.iter() { + let c = c as usize; + let (Some(m), Some(x), Some(mi)) = (main.faces.get(c), mx.get(c), info.get(c)) else { continue }; + if !mi.alive || m.corners().iter().any(|k| f.corners().contains(k)) { + continue; + } + if !(mi.length > 0.0) || dot(mi.area, normal) >= -FACING * mi.length { + continue; + } + let theirs: Vec = x.iter().map(|p| flat(*p)).collect(); + for s in &own_triangles { + for tm in triangles(m.count) { + let t = counter_clockwise(tm, &theirs); + let overlap = clip(&t, s); + if overlap.len() < 3 { + continue; + } + let Some(o) = overlap.first().copied() else { continue }; + for w in overlap.windows(2).skip(1) { + let (p, q) = (w[0], w[1]); + let tri_area = 0.5 * cross2(o, p, q); + if !(tri_area > 0.0) { + continue; + } + for (bary, weight) in TRIANGLE_RULE { + let y = [0, 1].map(|k| o[k] + bary[0] * (p[k] - o[k]) + bary[1] * (q[k] - o[k])); + let point = lift(y); + let (Some(qs), Some(qm)) = (closest_point(fx, point), closest_point(x, point)) else { continue }; + let gap = dot(sub(qs.point, qm.point), qm.normal); + // Over the main face (its foot not clamped to an edge), behind it, within + // its depth. + if !(gap < 0.0 && -gap <= m.depth && qm.lateral <= mi.edge) { + continue; + } + let push = eps * (-gap) * weight * tri_area; + if !push.is_finite() { + continue; + } + let force = qm.normal.map(|k| push * k); + for (node, w) in f.corners().iter().zip(qs.weights) { + out.forces.push((*node, force.map(|k| w * k))); + out.normal.push((*node, w * push, c)); + } + for (node, w) in m.corners().iter().zip(qm.weights) { + out.forces.push((*node, force.map(|k| -w * k))); + } + out.penetration = out.penetration.max(-gap); + out.points += 1; + out.total += push; + } + } + } + } + } + out +} + +impl Interface { + /// The contact `name` of the model on the analysed nodes and elements: its surfaces, and the + /// initial penetrations adjusted in `xyz` (or refused). + pub fn new( + model: &Model, + resolved: &Resolved, + name: &str, + contact: &Contact, + xyz: &mut [Vec3], + labels: &[u32], + elements: &[u32], + ) -> Result { + let context = |e: Error| match e { + Error::Invalid(m) => invalid(format!("contact `{name}`: {m}")), + other => other, + }; + let node = |l: u32| labels.binary_search(&l).ok(); + let element = |l: u32| elements.binary_search(&l).ok(); + let faces_of = |set: &str| resolved.sets.get(set).map(|s| s.faces.clone()).unwrap_or_default(); + let (surfaces, passes) = match &contact.surfaces { + ContactSurfaces::Pair { main, secondary } => { + let main = Surface::new(model, &faces_of(main), xyz, &node, &element).map_err(context)?; + let secondary = Surface::new(model, &faces_of(secondary), xyz, &node, &element).map_err(context)?; + (vec![main, secondary], vec![(1, 0), (0, 1)]) + } + ContactSurfaces::General { surface } => { + let faces = match surface { + Some(s) => faces_of(s), + None => model.exterior_faces(), + }; + (vec![Surface::new(model, &faces, xyz, &node, &element).map_err(context)?], vec![(0, 0)]) + } + }; + let passes = passes + .into_iter() + .map(|(secondary, main)| Pass { secondary, main, slip: vec![Slip::default(); surfaces.get(secondary).map_or(0, |s| s.nodes.len())] }) + .collect(); + let mut interface = Interface { + name: name.to_string(), + friction: contact.friction, + scale: contact.penalty_scale, + surfaces, + passes, + adjusted: Vec::new(), + mass_per_area: None, + }; + interface.adjust(xyz, labels, contact.initial_penetration)?; + Ok(interface) + } + + /// Move the nodes that start inside a face onto it, or refuse them. + fn adjust(&mut self, xyz: &mut [Vec3], labels: &[u32], mode: InitialPenetration) -> Result<()> { + let zero: Vec = Vec::new(); + for pass in &self.passes { + let (Some(nodes), Some(faces)) = (self.surfaces.get(pass.secondary), self.surfaces.get(pass.main)) else { continue }; + let positions = faces.positions(xyz, &zero); + let grid = grid(faces, &positions); + let own = nodes.positions(xyz, &zero); + let normals = nodes.node_normals(&own); + // Each move is decided on the positions before any move of this pass. + let mut moves: Vec<(usize, Vec3, Adjustment)> = Vec::new(); + for (i, &n) in nodes.nodes.iter().enumerate() { + let x = xyz.get(n).copied().unwrap_or_default(); + let normal = normals.get(i).copied().unwrap_or_default(); + let Some((c, q)) = inside(faces, &positions, &grid, n, x, normal) else { continue }; + let Some(face) = faces.faces.get(c) else { continue }; + // Coordinates are exact to a few units in the last place: a node a hair behind a + // face it was meshed on is touching, not inside. + let scale = x.iter().fold(1.0f64, |m, c| m.max(c.abs())); + if -q.gap <= 1e-9 * scale { + continue; + } + let label = labels.get(n).copied().unwrap_or(0); + let adjustment = Adjustment { node: label, element: face.element, face: face.face, depth: -q.gap }; + if mode == InitialPenetration::Refuse { + return Err(invalid(format!( + "contact `{}`: node {label} starts {:.4e} inside face {} of element {} (move it out, or let the contact adjust it: initial_penetration \"adjust\")", + self.name, -q.gap, face.face, face.element + ))); + } + moves.push((n, q.point, adjustment)); + } + for (n, point, a) in moves { + if let Some(x) = xyz.get_mut(n) { + *x = point; + } + self.adjusted.push(a); + } + } + Ok(()) + } + + /// The number of faces. + pub fn faces(&self) -> usize { + self.surfaces.iter().map(|s| s.faces.len()).sum() + } + + /// Add the contact forces at displacements `u` into `force` (per degree of freedom), with + /// `du` the displacement increment that led to `u` over `dt_prev`, the stable increment + /// `dt` and the lumped masses per degree of freedom `mass`; update the friction state. + #[allow(clippy::too_many_arguments)] + pub fn forces( + &mut self, + xyz: &[Vec3], + u: &[f64], + du: &[f64], + dt_prev: f64, + dt: f64, + mass: &[f64], + state: &State, + force: &mut [f64], + ) -> Evaluation { + let mut e = Evaluation::default(); + let node_mass = |n: usize| (0..3).map(|k| mass.get(NODE_DOFS * n + k).copied().unwrap_or(0.0)).fold(f64::INFINITY, f64::min); + let mass_per_area = *self.mass_per_area.get_or_insert_with(|| { + self.surfaces + .iter() + .flat_map(|s| s.nodes.iter().zip(&s.area)) + .filter(|(_, a)| **a > 0.0) + .map(|(n, a)| node_mass(*n) / a) + .fold(f64::INFINITY, f64::min) + }); + // Every contact runs both ways, at half the penalty each. + let eps = 0.5 * self.scale * GAMMA * mass_per_area / (dt * dt); + if !(eps > 0.0 && eps.is_finite()) { + for pass in &mut self.passes { + pass.slip.fill(Slip::default()); + } + return e; + } + let alive = |f: &Face| state.elements.get(f.index).is_none_or(|s| s.failed.is_none()); + let positions: Vec>> = self.surfaces.iter().map(|s| s.positions(xyz, u)).collect(); + let grids: Vec = self.surfaces.iter().zip(&positions).map(|(s, p)| grid(s, p)).collect(); + let friction = self.friction; + for pass in &mut self.passes { + let (Some(secondary), Some(main), Some(sx), Some(mx), Some(grid)) = ( + self.surfaces.get(pass.secondary), + self.surfaces.get(pass.main), + positions.get(pass.secondary), + positions.get(pass.main), + grids.get(pass.main), + ) else { + continue; + }; + let info = face_info(main, mx, &alive); + let blocks = fieldcraft_parallel::map_grained(secondary.faces.len().div_ceil(BLOCK), 1, |b| { + let mut found = Vec::new(); + (b * BLOCK..((b + 1) * BLOCK).min(secondary.faces.len())) + .map(|i| match (secondary.faces.get(i), sx.get(i)) { + (Some(f), Some(fx)) => integrate_face(f, fx, alive(f), main, mx, &info, grid, eps, &mut found), + _ => FaceForces::default(), + }) + .collect::>() + }); + // Per secondary node: its normal force, and the main face that gave it the most. + let mut normal = vec![0.0; secondary.nodes.len()]; + let mut partner: Vec> = vec![None; secondary.nodes.len()]; + for face in blocks.iter().flatten() { + for (node, f) in &face.forces { + add_at(force, *node, f); + } + for (node, share, c) in &face.normal { + if let Ok(k) = secondary.nodes.binary_search(node) + && let (Some(slot), Some(best)) = (normal.get_mut(k), partner.get_mut(k)) + { + *slot += share; + if best.is_none_or(|(_, s)| *share > s) { + *best = Some((*c, *share)); + } + } + } + e.normal_force += face.total; + e.penetration = e.penetration.max(face.penetration); + e.active += face.points; + } + if friction.is_frictionless() { + continue; + } + // Friction at the secondary nodes, in node order, against the point under the node + // of the main face that pressed it most. + for (k, (&n, slip)) in secondary.nodes.iter().zip(pass.slip.iter_mut()).enumerate() { + let fn_ = normal.get(k).copied().unwrap_or(0.0); + let kt = eps * secondary.area.get(k).copied().unwrap_or(0.0); + let hit = partner + .get(k) + .copied() + .flatten() + .filter(|_| fn_ > 0.0 && kt > 0.0) + .and_then(|(c, _)| mx.get(c).and_then(|x| closest_point(x, position(xyz, u, n))).map(|q| (c, q))); + let Some((c, q)) = hit else { + *slip = Slip::default(); + continue; + }; + let Some(face) = main.faces.get(c) else { continue }; + let mut rel = shift(du, n); + for (m, w) in face.corners().iter().zip(q.weights) { + let d = shift(du, *m); + for c in 0..3 { + rel[c] -= w * d[c]; + } + } + let normal = q.normal; + let rn = dot(rel, normal); + let slide = [0, 1, 2].map(|c| rel[c] - rn * normal[c]); + let old = slip.force; + let on = dot(old, normal); + let trial: Vec3 = std::array::from_fn(|c| old[c] - on * normal[c] - kt * slide[c]); + let size = norm(trial); + let rate = if dt_prev > 0.0 { norm(slide) / dt_prev } else { 0.0 }; + let limit = friction.coefficient(rate) * fn_; + let tangential = if size <= limit { + trial + } else { + e.dissipation += limit * (size - limit) / kt; + trial.map(|c| c * limit / size) + }; + slip.force = tangential; + add_at(force, n, &tangential); + for (m, w) in face.corners().iter().zip(q.weights) { + add_at(force, *m, &tangential.map(|c| -w * c)); + } + } + } + e + } +} + +/// Every contact of the model, built on the analysed nodes and elements (with the initial +/// adjustments made in `xyz`), and no more than `max_faces` faces in all. +pub fn interfaces( + model: &Model, + resolved: &Resolved, + xyz: &mut [Vec3], + labels: &[u32], + elements: &[u32], + max_faces: usize, +) -> Result> { + let mut out = Vec::with_capacity(model.contacts().len()); + let mut faces = 0usize; + for (name, c) in model.contacts() { + let i = Interface::new(model, resolved, name, c, xyz, labels, elements)?; + faces = faces.saturating_add(i.faces()); + if faces > max_faces { + return Err(Error::Limit(format!( + "the contacts hold more than {max_faces} faces (contact `{name}` passes the cap); name smaller surfaces" + ))); + } + out.push(i); + } + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn clipping_two_triangles_gives_their_overlap() { + let a = [[0.0, 0.0], [2.0, 0.0], [0.0, 2.0]]; + let b = [[1.0, -1.0], [1.0, 3.0], [-1.0, 1.0]]; + // b clockwise: reversed it is counter-clockwise. + let mut b2 = b.to_vec(); + if signed_area(&b2) < 0.0 { + b2.reverse(); + } + let o = clip(&a, &b2); + // The overlap: the part of a with x ≤ 1 and y ≤ x + 2… here the triangle's left half, + // area 2 − ½ = 1.5. + assert!((signed_area(&o) - 1.5).abs() < 1e-12, "{o:?} {}", signed_area(&o)); + assert!(clip(&a, &[[5.0, 5.0], [6.0, 5.0], [5.0, 6.0]]).is_empty()); + let same = clip(&a, &a); + assert!((signed_area(&same) - 2.0).abs() < 1e-12); + } +} diff --git a/crates/explicit/src/elements.rs b/crates/explicit/src/elements.rs index 3bbd6b4..17752d2 100644 --- a/crates/explicit/src/elements.rs +++ b/crates/explicit/src/elements.rs @@ -70,6 +70,8 @@ pub struct Workspace { pub stress_force: Vec, pub hourglass_force: Vec, pub viscous_force: Vec, + /// Contact forces on the model per degree of freedom (empty without contact). + pub contact_force: Vec, /// Per node: the volume-weighted volumetric strain increments of its tetrahedra, and their /// volume. theta: Vec, @@ -86,6 +88,7 @@ impl Workspace { stress_force: vec![0.0; dofs], hourglass_force: vec![0.0; dofs], viscous_force: vec![0.0; dofs], + contact_force: if sys.contacts.is_empty() { Vec::new() } else { vec![0.0; dofs] }, theta: vec![0.0; n], tet_volume: vec![0.0; n], tet_increment: vec![None; sys.tets.len()], diff --git a/crates/explicit/src/integrate.rs b/crates/explicit/src/integrate.rs index b162626..e8e48fc 100644 --- a/crates/explicit/src/integrate.rs +++ b/crates/explicit/src/integrate.rs @@ -6,13 +6,16 @@ //! 1. u^n = u^{n−1} + Δt^{n−½} v^{n−½}; t^n = t^{n−1} + Δt^{n−½}; //! 2. internal forces f^n from u^n (stresses, hourglass stiffness) and v^{n−½} (hourglass and //! bulk viscosity), and the next increment Δt^{n+½} = scale · min_e Δt_e; -//! 3. v^{n+½} = v^{n−½} + Δt^n M⁻¹(f_ext − f^n), then the constraints: rigid walls set the normal -//! velocity of a node that would cross its wall so it lands on it (a frictionless kinematic -//! contact condition: T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier, W. -//! Kanoknukulchai, *A finite element method for a class of contact-impact problems*, CMAME -//! 8(3):249–276, 1976), and prescribed degrees of freedom take their velocity. The change a -//! constraint makes is the impulse of its reaction, r^n = M Δv / Δt^n; -//! 4. energies: the works of the internal and external forces by the trapezoidal rule, +//! 3. the contact forces c^n at u^n ([`crate::contact`], with the stable increment Δt^{n+½}); +//! v^{n+½} = v^{n−½} + Δt^n M⁻¹(f_ext + c^n − f^n), then the constraints: rigid walls set the +//! normal velocity of a node that would cross its wall so it lands on it (a kinematic contact +//! condition: T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier, W. Kanoknukulchai, +//! *A finite element method for a class of contact-impact problems*, CMAME 8(3):249–276, +//! 1976), with Coulomb friction as a tangential impulse at most μ times the normal one (the +//! tangential velocity stops, or loses μ Δv_n when it slides), and prescribed degrees of +//! freedom take their velocity. The change a constraint makes is the impulse of its reaction, +//! r^n = M Δv / Δt^n; +//! 4. energies: the works of the internal, contact and external forces by the trapezoidal rule, //! W += ½ Δu · (f^{n−1} + f^n) (exactly the strain energy for a linear elastic material), //! the constraints' work Δt^n r^n · ½(v^{n−½} + v^{n+½}) (exactly the kinetic energy their //! impulses add or remove), and the kinetic energy of v^n = v^{n−½} + ½Δt^{n−½} a^n. @@ -22,8 +25,9 @@ use fieldcraft_model::{ExplicitDynamic, Model, OutputVariable, Procedure, Step}; +use crate::contact::{Evaluation, Interface}; use crate::elements::{Motion, State, Workspace, internal_forces}; -use crate::result::{ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History}; +use crate::result::{ContactSummary, ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History}; use crate::system::{NODE_DOFS, System}; use crate::{Error, Limits, Result}; @@ -70,8 +74,10 @@ struct Output { } /// Velocity constraints over an increment of length `dt` (rigid walls when dt > 0, then -/// prescribed degrees of freedom). -fn constrain(sys: &System, u: &[f64], v: &mut [f64], dt: f64) { +/// prescribed degrees of freedom); returns the kinetic energy the walls' friction took (`mass` +/// per degree of freedom). +fn constrain(sys: &System, u: &[f64], v: &mut [f64], dt: f64, mass: &[f64]) -> f64 { + let mut dissipated = 0.0; if dt > 0.0 { for w in &sys.walls { for &node in &w.nodes { @@ -86,6 +92,9 @@ fn constrain(sys: &System, u: &[f64], v: &mut [f64], dt: f64) { *c += change * w.normal[k]; } } + if !w.friction.is_frictionless() { + dissipated += wall_friction(v, node, w, change, mass.get(at(0)).copied().unwrap_or(0.0)); + } } } } @@ -95,14 +104,75 @@ fn constrain(sys: &System, u: &[f64], v: &mut [f64], dt: f64) { *c = *p; } } + dissipated } -/// v_out = v_in + dt M⁻¹ (f_ext − f_int). +/// Coulomb friction of a node the wall `w` stopped with a normal velocity change `change`: its +/// tangential velocity changes by at most μ · change (stick when that stops it); returns the +/// kinetic energy taken, ½ m (|v_t|² − |v_t′|²). +fn wall_friction(v: &mut [f64], node: usize, w: &crate::system::Wall, change: f64, mass: f64) -> f64 { + let at = |k: usize| NODE_DOFS * node + k; + let vel = [0, 1, 2].map(|k| v.get(at(k)).copied().unwrap_or(0.0)); + let vn: f64 = (0..3).map(|k| vel[k] * w.normal[k]).sum(); + let tangential = [0, 1, 2].map(|k| vel[k] - vn * w.normal[k]); + let slide = tangential.iter().map(|c| c * c).sum::().sqrt(); + if slide == 0.0 { + return 0.0; + } + let cap = w.friction.coefficient(slide) * change; + let keep = if slide <= cap { 0.0 } else { 1.0 - cap / slide }; + for k in 0..3 { + if let Some(c) = v.get_mut(at(k)) { + *c = vn * w.normal[k] + keep * tangential[k]; + } + } + 0.5 * mass * slide * slide * (1.0 - keep * keep) +} + +/// v_out = v_in + dt M⁻¹ (f_ext + f_contact − f_int). fn accelerate(sys: &System, ws: &Workspace, mass: &[f64], v_in: &[f64], dt: f64, v_out: &mut [f64]) { let internal = ws.stress_force.iter().zip(&ws.hourglass_force).zip(&ws.viscous_force); for ((((out, v), m), e), ((s, h), q)) in v_out.iter_mut().zip(v_in).zip(mass).zip(&sys.external).zip(internal) { *out = v + dt * ((e - (s + h + q)) / m); } + // Contact adds its forces (a separate pass, so a model without contact runs the same sums). + if !ws.contact_force.is_empty() { + for ((out, m), c) in v_out.iter_mut().zip(mass).zip(&ws.contact_force) { + *out += dt * (c / m); + } + } +} + +/// The contacts' forces at `u` (with the increment `du` over `dt_prev` that led there, and the +/// stable increment `dt`) into `ws.contact_force`, and what each found. +#[allow(clippy::too_many_arguments)] +fn contact_forces( + contacts: &mut [Interface], + sys: &System, + u: &[f64], + du: &[f64], + dt_prev: f64, + dt: f64, + mass: &[f64], + state: &State, + ws: &mut Workspace, +) -> Vec { + ws.contact_force.fill(0.0); + contacts.iter_mut().map(|c| c.forces(&sys.xyz, u, du, dt_prev, dt, mass, state, &mut ws.contact_force)).collect() +} + +/// A contact's record over the run. +fn track(summary: &mut ContactSummary, e: &Evaluation, t: f64) { + if e.active > 0 { + summary.first_contact.get_or_insert(t); + summary.last_contact = Some(t); + } + if e.normal_force > summary.peak_force { + summary.peak_force = e.normal_force; + summary.peak_time = t; + } + summary.max_penetration = summary.max_penetration.max(e.penetration); + summary.friction += e.dissipation; } /// The constraint reactions M (v_constrained − v_free) / Δt, into `r`. @@ -120,7 +190,7 @@ fn per_node(x: &[f64]) -> Vec<[f64; 3]> { pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Result { let resolved = model.resolve()?; let (step, settings) = select_step(model, step)?; - let sys = System::new(model, &resolved, step, settings, limits.max_dofs)?; + let mut sys = System::new(model, &resolved, step, settings, limits)?; let end = settings.end_time; let scale = settings.dt_scale; let bulk = settings.bulk_viscosity; @@ -192,6 +262,9 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res let total_mass: f64 = sys.mass.iter().sum::() + sys.discrete.mass.iter().sum::() / 3.0; let mut ws = Workspace::new(&sys); let mut previous = Workspace::new(&sys); + let mut contacts = std::mem::take(&mut sys.contacts); + let mut summaries: Vec = + contacts.iter().map(|c| ContactSummary { name: c.name.clone(), adjusted: c.adjusted.clone(), ..ContactSummary::default() }).collect(); let mut state = State::new(&sys); let mut eroded: Vec = Vec::new(); let mut u = vec![0.0; n]; @@ -205,10 +278,14 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res let crit = internal_forces(&sys, &start, &mut state, &mut ws)?; let mut dt_half = (scale * crit).min(end); let initial_increment = dt_half; + let mut found = contact_forces(&mut contacts, &sys, &u, &du, 0.0, scale * crit, &mass, &state, &mut ws); + for (summary, e) in summaries.iter_mut().zip(&found) { + track(summary, e, 0.0); + } let mut v_half = vec![0.0; n]; accelerate(&sys, &ws, &mass, &v0, 0.5 * dt_half, &mut v_half); let mut free = v_half.clone(); - constrain(&sys, &u, &mut v_half, dt_half); + let wall_friction = constrain(&sys, &u, &mut v_half, dt_half, &mass); let mut reaction = vec![0.0; n]; impulses(&mass, &free, &v_half, 0.5 * dt_half, &mut reaction); // The constraints' work is that of their impulses on the mean of the velocities before @@ -217,12 +294,18 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res let impulse_work = |r: &[f64], before: &[f64], after: &[f64], dt: f64| -> f64 { dt * r.iter().zip(before.iter().zip(after)).map(|(r, (a, b))| r * 0.5 * (a + b)).sum::() }; - let mut energy = Energies { kinetic: initial_kinetic, constraint: impulse_work(&reaction, &v0, &v_half, 0.5 * dt_half), ..Energies::default() }; + let mut energy = Energies { + kinetic: initial_kinetic, + constraint: impulse_work(&reaction, &v0, &v_half, 0.5 * dt_half), + friction: wall_friction, + ..Energies::default() + }; let mut history = History { nodes: sys.history_nodes.iter().filter_map(|i| sys.labels.get(*i).copied()).collect(), ..History::default() }; let mut out_frames = Vec::new(); - let sample = |history: &mut History, t: f64, dt: f64, e: Energies, u: &[f64], v: &[f64], r: &[f64]| { + let sample = |history: &mut History, t: f64, dt: f64, e: Energies, u: &[f64], v: &[f64], r: &[f64], found: &[Evaluation]| { history.time.push(t); + history.contact_force.push(found.iter().map(|e| e.normal_force).collect()); history.increment.push(dt); history.energies.push(e); let sum = |x: &[f64], weight: &dyn Fn(usize, usize) -> f64| { @@ -260,7 +343,7 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res let plastic = |state: &State| state.elements.iter().map(|e| e.plastic_work).sum::(); let mut histories = Schedule { interval: settings.history_interval(), next: settings.history_interval() }; let mut frames_due = Schedule { interval: settings.output_interval(), next: settings.output_interval() }; - sample(&mut history, 0.0, dt_half, energy, &u, &v0, &reaction); + sample(&mut history, 0.0, dt_half, energy, &u, &v0, &reaction, &found); out_frames.push(frame(0.0, &u, &v0, &state)); let mut t = 0.0; @@ -294,9 +377,16 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res let last = remaining <= 1e-12 * end; let dt_next = if last { 0.0 } else { (scale * crit).min(remaining) }; let dt_n = 0.5 * (dt_half + dt_next); + if !contacts.is_empty() { + found = contact_forces(&mut contacts, &sys, &u, &du, dt_half, scale * crit, &mass, &state, &mut ws); + for (summary, e) in summaries.iter_mut().zip(&found) { + track(summary, e, t); + energy.friction += e.dissipation; + } + } accelerate(&sys, &ws, &mass, &v_half, dt_n, &mut v_new); free.clone_from(&v_new); - constrain(&sys, &u, &mut v_new, dt_next); + energy.friction += constrain(&sys, &u, &mut v_new, dt_next, &mass); impulses(&mass, &free, &v_new, dt_n, &mut reaction); // One pass over the degrees of freedom: v^n = v^{n−½} + ½Δt^{n−½} a^n with the constrained // a^n = (v^{n+½} − v^{n−½}) / Δt^n, and the works (trapezoidal over the displacement for @@ -321,15 +411,22 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res energy.external += external; energy.constraint += dt_n * constraint; energy.kinetic = 0.5 * kinetic; + if !ws.contact_force.is_empty() { + let work: f64 = du.iter().zip(ws.contact_force.iter().zip(&previous.contact_force)).map(|(d, (c, c0))| d * (c0 + c)).sum(); + energy.contact -= 0.5 * work; + } peak_internal = peak_internal.max(energy.internal.abs()); - energy.balance = energy.kinetic + energy.internal + energy.hourglass + energy.viscous - energy.external - energy.constraint - initial_kinetic; + energy.balance = energy.kinetic + energy.internal + energy.hourglass + energy.viscous + energy.contact + - energy.external + - energy.constraint + - initial_kinetic; // Every element eroded leaves an infinite critical increment: the nodes fly freely. if !(energy.balance.is_finite() && !crit.is_nan()) { return Err(Error::Unstable(format!("the solution is not finite at time {t:.6e} (increment {increments})"))); } if histories.due(t) || last { energy.plastic = plastic(&state); - sample(&mut history, t, dt_half, energy, &u, &v_now, &reaction); + sample(&mut history, t, dt_half, energy, &u, &v_now, &reaction, &found); } if frames_due.due(t) || last { out_frames.push(frame(t, &u, &v_now, &state)); @@ -342,7 +439,7 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res } if !completed { if history.time.last() != Some(&t) { - sample(&mut history, t, dt_half, energy, &u, &v_now, &reaction); + sample(&mut history, t, dt_half, energy, &u, &v_now, &reaction, &found); } if out_frames.last().map(|f| f.time) != Some(t) { out_frames.push(frame(t, &u, &v_now, &state)); @@ -385,6 +482,7 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res energy, initial_kinetic, eroded, + contacts: summaries, warnings, }) } diff --git a/crates/explicit/src/lib.rs b/crates/explicit/src/lib.rs index c9aa5da..42f0cb0 100644 --- a/crates/explicit/src/lib.rs +++ b/crates/explicit/src/lib.rs @@ -14,8 +14,12 @@ //! hourglass control (viscous or stiffness), the Bonet–Burton average nodal pressure //! tetrahedron, von Neumann–Richtmyer/Wilkins bulk viscosity, and the element critical time //! increments. -//! - **Integration** ([`integrate`]): the time loop, kinematic constraints, energy accounting -//! with a balance check, history samples and field frames ([`result`]). +//! - **Contact** ([`contact`]): penalty contact between element faces (pairs both ways, single +//! surface and general contact) with Coulomb friction by a regularised return map, a bucket +//! grid search, and the stiffness from the masses and the stable increment. +//! - **Integration** ([`integrate`]): the time loop, kinematic constraints (rigid walls with +//! friction), energy accounting with a balance check, history samples and field frames +//! ([`result`]). //! //! With nonlinear geometry (the default) the elements are updated Lagrangian and the material //! state is integrated in rate form with an objective stress rate (Jaumann by Hughes–Winget, @@ -33,6 +37,7 @@ // Small fixed-size vector kernels read as the textbook formulas with index loops. #![allow(clippy::needless_range_loop)] +pub mod contact; pub mod discrete; pub mod elements; pub mod integrate; @@ -41,7 +46,7 @@ pub mod system; pub use elements::Failure; pub use integrate::solve_explicit; -pub use result::{ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History}; +pub use result::{ContactSummary, ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History}; #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum Error { @@ -82,16 +87,24 @@ pub struct Limits { pub max_increments: u64, /// Numbers kept in field frames and histories together (8 bytes each). pub max_output_values: usize, + /// Element faces of every contact surface together. + pub max_contact_faces: usize, } impl Limits { pub const MAX_DOFS: usize = 6_000_000; pub const MAX_INCREMENTS: u64 = fieldcraft_model::steps::MAX_INCREMENTS; pub const MAX_OUTPUT_VALUES: usize = 50_000_000; + pub const MAX_CONTACT_FACES: usize = 10_000_000; } impl Default for Limits { fn default() -> Limits { - Limits { max_dofs: Limits::MAX_DOFS, max_increments: Limits::MAX_INCREMENTS, max_output_values: Limits::MAX_OUTPUT_VALUES } + Limits { + max_dofs: Limits::MAX_DOFS, + max_increments: Limits::MAX_INCREMENTS, + max_output_values: Limits::MAX_OUTPUT_VALUES, + max_contact_faces: Limits::MAX_CONTACT_FACES, + } } } diff --git a/crates/explicit/src/result.rs b/crates/explicit/src/result.rs index 70bc4b1..9067a43 100644 --- a/crates/explicit/src/result.rs +++ b/crates/explicit/src/result.rs @@ -22,21 +22,34 @@ pub struct Energies { pub external: f64, /// Work of the boundary conditions and rigid walls on the model. pub constraint: f64, + /// Energy taken by contact: the work the model does against the contact forces (the penalty + /// energy stored in the contact springs, plus what friction dissipates there). + pub contact: f64, + /// Friction dissipation, in contact and on rigid walls (part of `contact` and of + /// `constraint`). + pub friction: f64, /// Plastic work: the part of `internal` dissipated by plastic flow. pub plastic: f64, /// Internal energy of the eroded elements when they eroded (part of `internal`). pub eroded: f64, - /// kinetic + internal + hourglass + viscous − external − constraint − the initial kinetic - /// energy: zero when energy is conserved. + /// kinetic + internal + hourglass + viscous + contact − external − constraint − the initial + /// kinetic energy: zero when energy is conserved. pub balance: f64, } impl Energies { /// The balance relative to the largest energy in the account. pub fn relative_error(&self, initial_kinetic: f64) -> f64 { - let scale = [initial_kinetic, self.kinetic, self.internal + self.hourglass + self.viscous, self.external.abs(), self.constraint.abs()] - .into_iter() - .fold(0.0, f64::max); + let scale = [ + initial_kinetic, + self.kinetic, + self.internal + self.hourglass + self.viscous, + self.external.abs(), + self.constraint.abs(), + self.contact.abs(), + ] + .into_iter() + .fold(0.0, f64::max); if scale > 0.0 { self.balance.abs() / scale } else { 0.0 } } } @@ -84,12 +97,33 @@ pub struct History { pub momentum: Vec<[f64; 3]>, /// Sum of the forces the boundary conditions and rigid walls exert on the model. pub constraint_force: Vec<[f64; 3]>, + /// Per contact (in [`ExplicitResult::contacts`] order): the normal force across it (the sum + /// of its constraints' normal forces). + pub contact_force: Vec>, /// The nodes with histories, and per sample their displacement and velocity (in this order). pub nodes: Vec, pub displacement: Vec>, pub velocity: Vec>, } +/// What a contact did over the run. +#[derive(Debug, Clone, Default, PartialEq, Serialize)] +pub struct ContactSummary { + pub name: String, + /// Nodes moved onto the other surface before the analysis. + pub adjusted: Vec, + /// The largest normal force across it, and when. + pub peak_force: f64, + pub peak_time: f64, + /// The first and the last time it carried a force (none when it never did). + pub first_contact: Option, + pub last_contact: Option, + /// The deepest penetration of a node into a face. + pub max_penetration: f64, + /// Friction dissipated in it. + pub friction: f64, +} + /// An element that limits the stable time increment. #[derive(Debug, Clone, Copy, PartialEq, Serialize)] pub struct ControllingElement { @@ -122,5 +156,7 @@ pub struct ExplicitResult { pub initial_kinetic: f64, /// Elements eroded during the run, in the order they eroded. pub eroded: Vec, + /// Per contact of the model, in name order. + pub contacts: Vec, pub warnings: Vec, } diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs index 897eb47..b0e6405 100644 --- a/crates/explicit/src/system.rs +++ b/crates/explicit/src/system.rs @@ -8,12 +8,13 @@ use fieldcraft_fem::plasticity::J2; use fieldcraft_fem::uniform::{Hex8, Tet4, modulus_bound}; use fieldcraft_fem::{DofKind, DofLayout, DofSet, face}; use fieldcraft_model::{ - Elastic, ElementKind, ExplicitDynamic, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step, StressRate, Target, + Elastic, ElementKind, ExplicitDynamic, Friction, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step, StressRate, Target, }; +use crate::contact::Interface; use crate::discrete::{Discrete, Rod}; use crate::elements::{HexScales, add_at, frequency_bound, hex_scales}; -use crate::{Error, Result}; +use crate::{Error, Limits, Result}; /// Degrees of freedom per node: the three translations, node-major (see [`System`]). pub const NODE_DOFS: usize = 3; @@ -81,13 +82,14 @@ pub struct Tet { pub thermal_rate: f64, } -/// A fixed, frictionless rigid plane. +/// A fixed rigid plane. #[derive(Debug, Clone)] pub struct Wall { pub point: [f64; 3], /// Unit normal, towards the side the nodes stay on. pub normal: [f64; 3], pub nodes: Vec, + pub friction: Friction, } /// The discrete system of an explicit dynamic step. @@ -117,6 +119,8 @@ pub struct System { /// Prescribed velocity per degree of freedom (`None` = free). pub prescribed: Vec>, pub walls: Vec, + /// The model's contacts ([`crate::contact`]). + pub contacts: Vec, /// Velocity per degree of freedom at the start (prescribed velocities included). pub initial_velocity: Vec, /// Constant external force per degree of freedom. @@ -192,16 +196,12 @@ fn check_formulation(kind: ElementKind, section_name: &str, f: &fieldcraft_model } impl System { - pub fn new(model: &Model, resolved: &Resolved, step: &Step, settings: &ExplicitDynamic, max_dofs: usize) -> Result { + pub fn new(model: &Model, resolved: &Resolved, step: &Step, settings: &ExplicitDynamic, limits: &Limits) -> Result { + let max_dofs = limits.max_dofs; if let Some(name) = model.couplings().keys().next() { return Err(Error::NotSupported(format!("couplings (coupling `{name}`) in explicit dynamics"))); } - if let Some(name) = model.contacts().keys().next() { - return Err(Error::NotSupported(format!("contact (`{name}`) in explicit dynamics"))); - } - if let Some(w) = settings.rigid_walls.iter().find(|w| !w.friction.is_frictionless()) { - return Err(Error::NotSupported(format!("friction on rigid walls (the wall on set `{}`)", w.set))); - } + // Elements by label, with their section's material. let mut materials: Vec = Vec::new(); let mut chosen: Vec<(u32, ElementKind, Vec, usize, Option)> = Vec::new(); @@ -274,6 +274,7 @@ impl System { discrete_increment: None, prescribed: vec![None; dofs], walls: Vec::new(), + contacts: Vec::new(), initial_velocity: vec![0.0; dofs], external: vec![0.0; dofs], history_nodes: Vec::new(), @@ -285,6 +286,20 @@ impl System { labels, xyz, }; + // Contact first: nodes that start inside a surface are moved before the elements are + // built on the coordinates. + sys.contacts = crate::contact::interfaces(model, resolved, &mut sys.xyz, &sys.labels, &sys.elements, limits.max_contact_faces)?; + for c in &sys.contacts { + if let Some(deepest) = c.adjusted.iter().max_by(|a, b| a.depth.total_cmp(&b.depth)) { + sys.warnings.push(format!( + "contact `{}`: {} nodes started inside the other surface and were moved onto it (the deepest, node {}, by {:.3e})", + c.name, + c.adjusted.len(), + deepest.node, + deepest.depth + )); + } + } for (index, (label, kind, nodes, mi, area)) in chosen.into_iter().enumerate() { let idx: Vec = nodes.iter().map(|l| sys.node(*l).ok_or_else(|| Error::Internal(format!("node {l}")))).collect::>()?; let mat = sys.materials.get(mi).ok_or_else(|| Error::Internal("material index".into()))?; @@ -590,7 +605,7 @@ impl System { return Err(invalid(format!("node {label} of set `{}` starts behind its rigid wall (gap {gap:e})", w.set))); } } - self.walls.push(Wall { point: w.point, normal, nodes }); + self.walls.push(Wall { point: w.point, normal, nodes, friction: w.friction }); } if let Some(set) = &settings.history_nodes { let nodes = self.nodes_of(resolved, set); diff --git a/crates/explicit/tests/contact.rs b/crates/explicit/tests/contact.rs new file mode 100644 index 0000000..822a7eb --- /dev/null +++ b/crates/explicit/tests/contact.rs @@ -0,0 +1,707 @@ +//! Deformable contact and friction in explicit dynamics against closed-form solutions: +//! - two identical elastic bars colliding exchange their velocities after a contact of 2L/c, with +//! the interface force ρ c A v₀/2 (one-dimensional wave theory: K. F. Graff, *Wave Motion in +//! Elastic Solids*, Dover 1991, §2.1–2.3; W. Johnson, *Impact Strength of Materials*, Arnold +//! 1972, §1.6), and keep their energy (frictionless elastic impact); +//! - a block sliding on a plane with Coulomb friction decelerates at μg and stops after +//! v₀²/(2μg), on a rigid wall and on a deformable plate (any dynamics text, e.g. J. L. Meriam, +//! L. G. Kraige, *Engineering Mechanics: Dynamics*, §3/5), and with the exponential decay law +//! after the distance the decaying coefficient gives; +//! - two elastic spheres colliding: Hertz's contact duration and peak force (H. Hertz, *Über die +//! Berührung fester elastischer Körper*, J. reine angew. Math. 92:156–171, 1882; K. L. Johnson, +//! *Contact Mechanics*, Cambridge 1985, §11.4); +//! - a contact patch test: two blocks with non-matching meshes pressed together carry a uniform +//! stress across the interface (R. L. Taylor, P. Papadopoulos, *On a patch test for contact +//! problems in two dimensions*, in Nonlinear Computational Mechanics, Springer 1991, 690–702; +//! T. Y. Crisfield-style uniform-pressure form for three dimensions); +//! - and what a contact must never do: hold back a body in free flight, depend on the thread +//! count, or crash on hostile input. + +use fieldcraft_explicit::{ExplicitResult, Limits, solve_explicit}; +use fieldcraft_model::*; + +/// A model built block by block, each block of hex8 with its own nodes. +struct Mesh { + model: Model, + next_node: u32, + next_element: u32, +} + +impl Mesh { + fn new() -> Mesh { + Mesh { model: Model::default(), next_node: 1, next_element: 1 } + } + + /// A box from `origin` of `size` with `n` hex8 along each axis, positions mapped by `map`; + /// its elements. + fn block_mapped(&mut self, origin: [f64; 3], size: [f64; 3], n: [usize; 3], map: impl Fn([f64; 3]) -> [f64; 3]) -> Vec { + let base = self.next_node; + let at = |i: usize, j: usize, k: usize| base + (i + (n[0] + 1) * (j + (n[1] + 1) * k)) as u32; + for k in 0..=n[2] { + for j in 0..=n[1] { + for i in 0..=n[0] { + let x = [0, 1, 2].map(|d| origin[d] + size[d] * [i, j, k][d] as f64 / n[d] as f64); + self.model.set_node(at(i, j, k), map(x)).unwrap(); + } + } + } + self.next_node = at(n[0], n[1], n[2]) + 1; + let mut out = Vec::new(); + for k in 0..n[2] { + for j in 0..n[1] { + for i in 0..n[0] { + let nodes = vec![ + at(i, j, k), + at(i + 1, j, k), + at(i + 1, j + 1, k), + at(i, j + 1, k), + at(i, j, k + 1), + at(i + 1, j, k + 1), + at(i + 1, j + 1, k + 1), + at(i, j + 1, k + 1), + ]; + self.model.set_element(self.next_element, ElementKind::Hex8, nodes).unwrap(); + out.push(self.next_element); + self.next_element += 1; + } + } + } + out + } + + fn block(&mut self, origin: [f64; 3], size: [f64; 3], n: [usize; 3]) -> Vec { + self.block_mapped(origin, size, n, |x| x) + } + + fn xyz(&self, node: u32) -> [f64; 3] { + self.model.node(node).unwrap().xyz + } + + /// The faces of `elements` whose centre satisfies `f`. + fn faces_where(&self, elements: &[u32], f: impl Fn([f64; 3]) -> bool) -> Vec { + let mut out = Vec::new(); + for &e in elements { + let el = self.model.element(e).unwrap(); + for (i, local) in el.kind.faces().iter().enumerate() { + let c = [0, 1, 2].map(|d| local.iter().map(|n| self.xyz(el.nodes[*n])[d]).sum::() / local.len() as f64); + if f(c) { + out.push(ElementFace { element: e, face: i as u8 + 1 }); + } + } + } + out + } + + fn nodes_where(&self, elements: &[u32], f: impl Fn([f64; 3]) -> bool) -> Vec { + let mut out: Vec = elements.iter().flat_map(|e| self.model.element(*e).unwrap().nodes.clone()).filter(|n| f(self.xyz(*n))).collect(); + out.sort_unstable(); + out.dedup(); + out + } + + fn material(&mut self, name: &str, young: f64, poisson: f64, density: f64) { + let m = Material { elastic: Some(Elastic::Isotropic { young, poisson }), density: Some(density), ..Material::default() }; + self.model.add_material(name, m).unwrap(); + } + + /// A solid section of one-point hex8 of `material` on `elements`, under the set `name`. + fn section(&mut self, name: &str, elements: &[u32], material: &str) { + self.model.add_set(name, Set { elements: elements.to_vec(), ..Set::default() }).unwrap(); + let formulation = Formulation { integration: Integration::Reduced, ..Formulation::default() }; + self.model.add_section(name, Section { set: name.into(), material: material.into(), kind: SectionKind::Solid {}, formulation }).unwrap(); + } + + fn face_set(&mut self, name: &str, faces: Vec) { + assert!(!faces.is_empty(), "{name}"); + self.model.add_set(name, Set { faces, ..Set::default() }).unwrap(); + } + + fn node_set(&mut self, name: &str, nodes: Vec) { + assert!(!nodes.is_empty(), "{name}"); + self.model.add_set(name, Set { nodes, ..Set::default() }).unwrap(); + } + + fn pair(&mut self, name: &str, main: &str, secondary: &str, friction: Friction) { + let c = Contact { friction, ..Contact::new(ContactSurfaces::Pair { main: main.into(), secondary: secondary.into() }) }; + self.model.add_contact(name, c).unwrap(); + } + + fn step(&mut self, step: ExplicitDynamic, bcs: &[(&str, BoundaryCondition)], loads: &[(&str, Load)]) { + let mut s = Step::new("Run", Procedure::ExplicitDynamic(step)); + for (name, bc) in bcs { + self.model.add_boundary_condition(name, bc.clone()).unwrap(); + s.boundary_conditions.insert(name.to_string()); + } + for (name, load) in loads { + self.model.add_load(name, load.clone()).unwrap(); + s.loads.insert(name.to_string()); + } + self.model.add_step(s, None).unwrap(); + } + + fn solve(&self) -> ExplicitResult { + let r = solve_explicit(&self.model, Some("Run"), &Limits::default()).unwrap(); + assert!(r.completed, "{:?}", r.warnings); + r + } +} + +fn node_index(r: &ExplicitResult, label: u32) -> usize { + r.nodes.binary_search(&label).unwrap() +} + +/// The mean velocity of some nodes in the last frame. +fn mean_velocity(r: &ExplicitResult, nodes: &[u32]) -> [f64; 3] { + let f = r.frames.last().unwrap(); + let mut v = [0.0; 3]; + for n in nodes { + let x = f.velocity[node_index(r, *n)]; + for k in 0..3 { + v[k] += x[k] / nodes.len() as f64; + } + } + v +} + +const BAR_E: f64 = 1.0e4; +const BAR_RHO: f64 = 1.0; +const BAR_L: f64 = 10.0; +const BAR_W: f64 = 0.1; + +/// Two bars of `cells` hex8, the first moving at v₀ along x into the second at rest (touching +/// end to end, separate nodes), with `contact` between their ends: (mesh, bar A, bar B). +fn bars(cells: usize, v0: f64, end: f64, contact: impl FnOnce(&mut Mesh, Vec, Vec)) -> (Mesh, Vec, Vec) { + let mut m = Mesh::new(); + let a = m.block([0.0; 3], [BAR_L, BAR_W, BAR_W], [cells, 1, 1]); + let b = m.block([BAR_L, 0.0, 0.0], [BAR_L, BAR_W, BAR_W], [cells, 1, 1]); + m.material("Mat", BAR_E, 0.0, BAR_RHO); + m.section("A", &a, "Mat"); + m.section("B", &b, "Mat"); + let end_a = m.faces_where(&a, |c| (c[0] - BAR_L).abs() < 1e-9); + let start_b = m.faces_where(&b, |c| (c[0] - BAR_L).abs() < 1e-9); + contact(&mut m, end_a, start_b); + let mut step = ExplicitDynamic::new(end); + step.initial_velocities = vec![InitialVelocity { set: "A".into(), velocity: [v0, 0.0, 0.0] }]; + // The elastic solution: no linear bulk viscosity (it dissipates at the wave fronts). + step.bulk_viscosity = BulkViscosity { linear: 0.0, quadratic: 1.5 }; + step.history_interval = Some(end / 4000.0); + m.step(step, &[], &[]); + (m, a, b) +} + +/// Velocity exchange, contact duration 2L/c, interface force ρ c A v₀ / 2 and the energy of a +/// frictionless elastic impact, for a pair (both ways) and for general contact. +#[test] +fn colliding_bars_exchange_their_velocities() { + let c = (BAR_E / BAR_RHO).sqrt(); + let v0 = 0.1; + let duration = 2.0 * BAR_L / c; + for kind in ["pair", "general"] { + let (m, a, b) = bars(200, v0, 1.5 * duration, |m, end_a, start_b| match kind { + "pair" => { + m.face_set("end_a", end_a); + m.face_set("start_b", start_b); + m.pair("impact", "start_b", "end_a", Friction::NONE); + } + _ => m.model.add_contact("everything", Contact::new(ContactSurfaces::General { surface: None })).unwrap(), + }); + let r = m.solve(); + let s = &r.contacts[0]; + let (first, last) = (s.first_contact.unwrap(), s.last_contact.unwrap()); + let plateau = BAR_RHO * c * BAR_W * BAR_W * v0 / 2.0; + let va = mean_velocity(&r, &m.nodes_where(&a, |_| true))[0]; + let vb = mean_velocity(&r, &m.nodes_where(&b, |_| true))[0]; + let ke = r.energy.kinetic / r.initial_kinetic; + // The force across the interface while the waves cross the bars. + let mid: Vec = r + .history + .time + .iter() + .zip(&r.history.contact_force) + .filter(|(t, _)| **t > 0.2 * duration && **t < 0.8 * duration) + .map(|(_, f)| f[0]) + .collect(); + let mean_force = mid.iter().sum::() / mid.len() as f64; + eprintln!( + "{kind}: contact {first:.5}..{last:.5} ({:.5}) vs 2L/c = {duration:.5}; force {mean_force:.5e} (peak {:.5e}) vs ρcAv₀/2 = {plateau:.5e}; v_A {va:.5} v_B {vb:.5} vs 0, {v0}; KE/KE₀ {ke:.5}; contact energy {:.3e} of KE₀ {:.3e}; balance {:.2e}; max penetration {:.3e}; {} increments", + last - first, + s.peak_force, + r.energy.contact, + r.initial_kinetic, + r.energy.relative_error(r.initial_kinetic), + s.max_penetration, + r.increments + ); + assert!(first < 0.01 * duration, "{kind}: contact from the start, not {first}"); + assert!(((last - first) - duration).abs() < 0.015 * duration, "{kind}: contact for {} vs {duration}", last - first); + assert!((mean_force - plateau).abs() < 0.02 * plateau, "{kind}: {mean_force} vs {plateau}"); + assert!(va.abs() < 0.02 * v0 && (vb - v0).abs() < 0.02 * v0, "{kind}: {va} {vb}"); + assert!(ke > 0.97 && ke <= 1.0 + 1e-6, "{kind}: KE ratio {ke}"); + // After the separation the contact springs are empty again: what they took they gave back. + assert!(r.energy.contact.abs() < 0.005 * r.initial_kinetic, "{kind}: {:?}", r.energy); + assert!(r.energy.relative_error(r.initial_kinetic) < 0.005, "{kind}: {:?}", r.energy); + assert!(s.max_penetration < 0.01 * BAR_L / 200.0 * 10.0, "{kind}: penetration {}", s.max_penetration); + assert_eq!(r.energy.friction, 0.0); + } +} + +/// A body in free flight with general contact on all its faces feels no contact force (no node +/// grabs a face of its own body), and its motion is that of the same run without contact. +#[test] +fn general_contact_leaves_a_free_body_alone() { + let run = |contact: bool| { + let mut m = Mesh::new(); + let e = m.block([0.0; 3], [1.0, 0.6, 0.4], [5, 3, 2]); + m.material("Mat", 1e4, 0.3, 1.0); + m.section("all", &e, "Mat"); + if contact { + m.model.add_contact("all", Contact::new(ContactSurfaces::General { surface: None })).unwrap(); + } + let mut step = ExplicitDynamic::new(0.05); + step.initial_velocities = vec![InitialVelocity { set: "all".into(), velocity: [1.0, -0.5, 0.2] }]; + m.step(step, &[], &[]); + m.solve() + }; + let (with, without) = (run(true), run(false)); + assert!(with.history.contact_force.iter().all(|f| f[0] == 0.0)); + assert_eq!(with.contacts[0].first_contact, None); + assert_eq!(with.frames.last().unwrap().displacement, without.frames.last().unwrap().displacement); + assert_eq!(with.energy.contact, 0.0); +} + +/// What kind of plane a block slides on. +#[derive(Clone, Copy, PartialEq, Debug)] +enum Plane { + RigidWall, + /// A plate fixed at its base, in contact with the block. + Plate, +} + +/// A block (1 × 1 × 0.5, 4 × 4 × 2 hex8) launched along x at v₀ on a horizontal plane under +/// gravity g with `friction`: the distance its centre of mass travels and the time it stops. +fn slide(plane: Plane, friction: Friction, v0: f64, g: f64, end: f64) -> (f64, f64, ExplicitResult) { + let mut m = Mesh::new(); + let block = m.block([0.0, 0.0, 0.0], [1.0, 1.0, 0.5], [4, 4, 2]); + m.material("Mat", 1.0e5, 0.3, 1.0); + m.section("block", &block, "Mat"); + let bottom = m.nodes_where(&block, |x| x[2].abs() < 1e-9); + let mut step = ExplicitDynamic::new(end); + step.initial_velocities = vec![InitialVelocity { set: "block".into(), velocity: [v0, 0.0, 0.0] }]; + step.history_interval = Some(end / 2000.0); + let mut bcs = Vec::new(); + match plane { + Plane::RigidWall => { + m.node_set("bottom", bottom); + step.rigid_walls = vec![RigidWall { set: "bottom".into(), point: [0.0; 3], normal: [0.0, 0.0, 1.0], friction }]; + } + Plane::Plate => { + let plate = m.block([-0.5, -0.5, -0.5], [4.5, 2.0, 0.5], [9, 4, 1]); + m.section("plate", &plate, "Mat"); + let top = m.faces_where(&plate, |c| c[2].abs() < 1e-9); + let under = m.faces_where(&block, |c| c[2].abs() < 1e-9); + m.face_set("top", top); + m.face_set("under", under); + m.pair("slide", "top", "under", friction); + m.node_set("base", m.nodes_where(&plate, |x| (x[2] + 0.5).abs() < 1e-9)); + bcs.push(("hold", BoundaryCondition::Fixed { target: Target::Set("base".into()) })); + } + } + let gravity = Load::Gravity { acceleration: [0.0, 0.0, -g], set: Some("block".into()) }; + m.step(step, &bcs, &[("weight", gravity)]); + let r = m.solve(); + let nodes = m.nodes_where(&block, |_| true); + let f = r.frames.last().unwrap(); + let distance = nodes.iter().map(|n| f.displacement[node_index(&r, *n)][0]).sum::() / nodes.len() as f64; + // The block's momentum per mass (the plate's is nothing along x when it is held). + let mass = match plane { + Plane::RigidWall => r.total_mass, + Plane::Plate => 0.5, + }; + let stop = r.history.time.iter().zip(&r.history.momentum).find(|(_, p)| p[0] / mass < 0.005 * v0).map_or(f64::NAN, |(t, _)| *t); + (distance, stop, r) +} + +/// Coulomb friction μ: deceleration μg, stopping after v₀/(μg) at v₀²/(2μg); on a rigid wall +/// (kinematic friction impulses) and on a deformable plate (penalty contact, return map). +#[test] +fn a_block_slides_to_a_stop_under_coulomb_friction() { + let (v0, g, mu) = (3.0, 10.0, 0.3); + let (t_stop, d_stop) = (v0 / (mu * g), v0 * v0 / (2.0 * mu * g)); + for plane in [Plane::RigidWall, Plane::Plate] { + let (d, t, r) = slide(plane, Friction::coulomb(mu), v0, g, 1.3 * t_stop); + eprintln!( + "{plane:?}: stops at {t:.4} vs v₀/(μg) = {t_stop:.4}, after {d:.4} vs v₀²/(2μg) = {d_stop:.4}; friction dissipated {:.4} vs ½mv₀² = {:.4}; balance {:.2e}; {} increments", + r.energy.friction, + 0.5 * 0.5 * v0 * v0, + r.energy.relative_error(r.initial_kinetic), + r.increments + ); + assert!((t - t_stop).abs() < 0.02 * t_stop, "{plane:?}: {t} vs {t_stop}"); + assert!((d - d_stop).abs() < 0.02 * d_stop, "{plane:?}: {d} vs {d_stop}"); + // The kinetic energy went into friction (and a little into the bouncing on the plane). + let ke0 = 0.5 * 0.5 * v0 * v0; + assert!((r.energy.friction - ke0).abs() < 0.03 * ke0, "{plane:?}: {:?}", r.energy); + assert!(r.energy.relative_error(r.initial_kinetic) < 0.01, "{plane:?}: {:?}", r.energy); + } +} + +/// The exponential decay law μ(v) = μ_k + (μ_s − μ_k) e^{−d v}: the block decelerates at +/// g μ(v), so it stops after t = ∫ dv / (g μ(v)) = [v₀ + ln((μ_k + Δμ e^{−d v₀}) / μ_s) / d] / (g μ_k) +/// and travels ∫ v dv / (g μ(v)) (by quadrature). +#[test] +fn a_block_slides_with_a_decaying_friction_coefficient() { + let (v0, g) = (3.0, 10.0); + let f = Friction { static_coefficient: 0.5, kinetic: Some(0.2), decay: 1.5 }; + let (a, b, d): (f64, f64, f64) = (0.2, 0.3, 1.5); + let t_stop = (v0 + ((a + b * (-d * v0).exp()) / (a + b)).ln() / d) / (g * a); + let n = 20_000; + let d_stop: f64 = (0..n) + .map(|i| { + let v = (i as f64 + 0.5) * v0 / n as f64; + v / (g * f.coefficient(v)) * v0 / n as f64 + }) + .sum(); + for plane in [Plane::RigidWall, Plane::Plate] { + let (dist, t, _) = slide(plane, f, v0, g, 1.3 * t_stop); + eprintln!("{plane:?}, decaying μ: stops at {t:.4} vs {t_stop:.4}, after {dist:.4} vs {d_stop:.4}"); + assert!((t - t_stop).abs() < 0.03 * t_stop, "{plane:?}: {t} vs {t_stop}"); + assert!((dist - d_stop).abs() < 0.03 * d_stop, "{plane:?}: {dist} vs {d_stop}"); + } +} + +/// Two blocks with non-matching meshes (3 × 3 below, 2 × 2 above) pressed together slowly +/// between symmetry planes: the stress is uniform σ_zz on both sides of the interface and the +/// force across it is σ_zz A (a contact patch test in uniaxial stress). +#[test] +fn a_contact_patch_test_carries_a_uniform_stress_across_non_matching_meshes() { + let mut m = Mesh::new(); + let lower = m.block([0.0; 3], [1.0, 1.0, 1.0], [3, 3, 2]); + let upper = m.block([0.0, 0.0, 1.0], [1.0, 1.0, 1.0], [2, 2, 2]); + let (young, rho) = (1000.0, 1.0); + m.material("Mat", young, 0.3, rho); + m.section("lower", &lower, "Mat"); + m.section("upper", &upper, "Mat"); + let both: Vec = lower.iter().chain(&upper).copied().collect(); + m.face_set("lower_top", m.faces_where(&lower, |c| (c[2] - 1.0).abs() < 1e-9)); + m.face_set("upper_bottom", m.faces_where(&upper, |c| (c[2] - 1.0).abs() < 1e-9)); + m.pair("interface", "lower_top", "upper_bottom", Friction::NONE); + m.node_set("x0", m.nodes_where(&both, |x| x[0].abs() < 1e-9)); + m.node_set("y0", m.nodes_where(&both, |x| x[1].abs() < 1e-9)); + m.node_set("z0", m.nodes_where(&both, |x| x[2].abs() < 1e-9)); + m.node_set("top", m.nodes_where(&both, |x| (x[2] - 2.0).abs() < 1e-9)); + let squeeze = 0.01; + let held = |set: &str, k: usize| { + let mut values = [None; 6]; + values[k] = Some(0.0); + BoundaryCondition::Displacement { target: Target::Set(set.into()), values } + }; + let mut push_values = [None; 6]; + push_values[2] = Some(-squeeze); + let push = BoundaryCondition::Displacement { target: Target::Set("top".into()), values: push_values }; + // A hundred wave transits: quasi-static. + let c = (young * 0.7 / (1.3 * 0.4) / rho).sqrt(); + let end = 100.0 * 2.0 / c; + let mut step = ExplicitDynamic::new(end); + step.output_interval = Some(end / 4.0); + m.step(step, &[("x0", held("x0", 0)), ("y0", held("y0", 1)), ("z0", held("z0", 2)), ("push", push)], &[]); + let r = m.solve(); + let f = r.frames.last().unwrap(); + let szz: Vec = f.stress.iter().map(|s| s[2]).collect(); + let mean = szz.iter().sum::() / szz.len() as f64; + let spread = szz.iter().map(|s| (s - mean).abs()).fold(0.0, f64::max); + let lateral = f.stress.iter().flat_map(|s| [s[0], s[1], s[3], s[4], s[5]]).map(f64::abs).fold(0.0, f64::max); + let force = *r.history.contact_force.last().unwrap().first().unwrap(); + eprintln!( + "patch: σ_zz mean {mean:.5} (E ε = {:.5}), spread {spread:.3e}, largest other component {lateral:.3e}; force across {force:.5} vs |σ_zz| A = {:.5}", + -young * squeeze / 2.0, + mean.abs() + ); + assert!(spread < 0.02 * mean.abs(), "σ_zz {szz:?}"); + assert!(lateral < 0.02 * mean.abs(), "{lateral}"); + assert!((force - mean.abs()).abs() < 0.02 * mean.abs(), "{force} vs {mean}"); + // The interface's penalty takes some of the squeeze (uniformly): each block is strained by + // what is left. + let closure = r.contacts[0].max_penetration; + let expected = young * (squeeze - closure) / 2.0; + eprintln!("patch: interface closure {closure:.4e}, σ_zz vs E (squeeze − closure) / 2 = {expected:.5}"); + assert!((mean.abs() - expected).abs() < 0.01 * expected, "{mean} vs {expected}"); +} + +/// The same contact run gives the same numbers, bit for bit, on one thread and on four. +#[test] +fn contact_results_do_not_depend_on_the_thread_count() { + let c = (BAR_E / BAR_RHO).sqrt(); + let (m, _, _) = bars(60, 0.1, 2.0 * BAR_L / c, |m, _, _| { + m.model.add_contact("all", Contact { friction: Friction::coulomb(0.2), ..Contact::new(ContactSurfaces::General { surface: None }) }).unwrap() + }); + let one = fieldcraft_parallel::with_threads(1, || solve_explicit(&m.model, None, &Limits::default()).unwrap()); + let four = fieldcraft_parallel::with_threads(4, || solve_explicit(&m.model, None, &Limits::default()).unwrap()); + assert_eq!(one.frames, four.frames); + assert_eq!(one.history, four.history); + assert_eq!(one.energy, four.energy); + assert_eq!(one.contacts, four.contacts); +} + +/// Nodes that start inside the other surface are moved onto it and reported, or refused with +/// the node, the face and the depth; a cap on contact faces is honoured; contact on the faces +/// of tetrahedra works too. +#[test] +fn initial_penetration_caps_and_tetrahedra() { + let build = |mode: InitialPenetration, overlap: f64| { + let mut m = Mesh::new(); + let a = m.block([0.0; 3], [1.0, 1.0, 1.0], [2, 2, 2]); + let b = m.block([0.0, 0.0, 1.0 - overlap], [1.0, 1.0, 1.0], [3, 3, 2]); + m.material("Mat", 1000.0, 0.3, 1.0); + m.section("a", &a, "Mat"); + m.section("b", &b, "Mat"); + m.face_set("a_top", m.faces_where(&a, |c| (c[2] - 1.0).abs() < 1e-9)); + m.face_set("b_bottom", m.faces_where(&b, |c| (c[2] - 1.0 + overlap).abs() < 1e-9)); + let c = Contact { initial_penetration: mode, ..Contact::new(ContactSurfaces::Pair { main: "a_top".into(), secondary: "b_bottom".into() }) }; + m.model.add_contact("stack", c).unwrap(); + m.step(ExplicitDynamic::new(0.01), &[], &[]); + m + }; + let r = build(InitialPenetration::Adjust, 0.01).solve(); + let s = &r.contacts[0]; + // Both passes: b's 16 bottom nodes into a, and a's 9 top nodes into b (they are inside b + // after b's nodes moved down? no: b's nodes are moved onto a's top first, then a's top nodes + // are on b's bottom). + assert!(!s.adjusted.is_empty(), "{s:?}"); + assert!(s.adjusted.iter().all(|a| (a.depth - 0.01).abs() < 1e-9), "{:?}", s.adjusted); + assert!(r.warnings.iter().any(|w| w.contains("contact `stack`") && w.contains("moved onto it")), "{:?}", r.warnings); + // After the adjustment nothing is pushed apart: no contact force worth the name. + assert!(s.peak_force < 1e-6, "{s:?}"); + let e = solve_explicit(&build(InitialPenetration::Refuse, 0.01).model, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("contact `stack`") && e.contains("starts 1.0000e-2 inside face 2 of element"), "{e}"); + // Touching exactly is not a penetration. + assert!(build(InitialPenetration::Refuse, 0.0).solve().contacts[0].adjusted.is_empty()); + let capped = Limits { max_contact_faces: 10, ..Limits::default() }; + let e = solve_explicit(&build(InitialPenetration::Adjust, 0.0).model, None, &capped).unwrap_err().to_string(); + assert!(e.contains("limit") && e.contains("contact `stack`"), "{e}"); +} + +/// Interval ends on [0, 1] in `n` intervals growing from the first by `ratio`. +fn graded(n: usize, ratio: f64) -> Vec { + let widths: Vec = (0..n).map(|k| ratio.powi(k as i32)).collect(); + let total: f64 = widths.iter().sum(); + let mut out = vec![0.0]; + let mut at = 0.0; + for w in widths { + at += w / total; + out.push(at); + } + if let Some(last) = out.last_mut() { + *last = 1.0; + } + out +} + +/// A quarter (x ≥ 0, y ≥ 0) of a ball of radius `radius` about `center`, meshed as a cubed sphere +/// of hex8: an inner cube of half-size 0.4 R with `n` uniform elements per edge, and six blocks +/// from its faces to the sphere with `n` × `n` elements across (equiangular, clustered towards +/// each face's centre by `lateral`) and `layers` elements through (thinner towards the surface +/// by `radial`). Its elements and the faces on the sphere within 30° of the pole below (−z) or +/// above (+z) the centre (`towards`). +fn quarter_ball( + m: &mut Mesh, + center: [f64; 3], + radius: f64, + n: usize, + layers: usize, + lateral: f64, + radial: f64, + towards: f64, +) -> (Vec, Vec) { + use std::collections::BTreeMap; + assert!(n.is_multiple_of(2)); + let inner = 0.4 * radius; + // Clustered lateral parameter on [−1, 1]: symmetric, finest at 0. + let half = graded(n / 2, lateral); + let mut clustered: Vec = half.iter().rev().map(|s| -s).collect(); + clustered.extend(half.iter().skip(1)); + let uniform: Vec = (0..=n).map(|i| -1.0 + 2.0 * i as f64 / n as f64).collect(); + // Radial blend from the cube (0) to the sphere (1), finest at the sphere. + let through: Vec = graded(layers, radial).iter().rev().map(|s| 1.0 - s).collect(); + let mut labels: BTreeMap<[i64; 3], u32> = BTreeMap::new(); + let quantum = 1e-9 * radius; + let mut node = |m: &mut Mesh, x: [f64; 3]| -> u32 { + let key = x.map(|c| (c / quantum).round() as i64); + *labels.entry(key).or_insert_with(|| { + let l = m.next_node; + m.model.set_node(l, [center[0] + x[0], center[1] + x[1], center[2] + x[2]]).unwrap(); + m.next_node += 1; + l + }) + }; + let kept = |x: &[[f64; 3]]| { + let c = [0, 1].map(|d| x.iter().map(|p| p[d]).sum::() / x.len() as f64); + c[0] > 0.0 && c[1] > 0.0 + }; + let mut hexes: Vec<[[f64; 3]; 8]> = Vec::new(); + // The inner cube. + let cube = |i: usize, j: usize, k: usize| [uniform[i], uniform[j], uniform[k]].map(|s| inner * s); + for k in 0..n { + for j in 0..n { + for i in 0..n { + hexes.push([ + cube(i, j, k), + cube(i + 1, j, k), + cube(i + 1, j + 1, k), + cube(i, j + 1, k), + cube(i, j, k + 1), + cube(i + 1, j, k + 1), + cube(i + 1, j + 1, k + 1), + cube(i, j + 1, k + 1), + ]); + } + } + } + // The six blocks: face axis a, sign σ; (u, v) along the other two axes. + for axis in 0..3 { + for sign in [-1.0, 1.0] { + let place = |u: f64, v: f64, w: f64| -> [f64; 3] { + let mut p = [0.0; 3]; + p[axis] = w; + p[(axis + 1) % 3] = u; + p[(axis + 2) % 3] = v; + p + }; + let point = |i: usize, j: usize, k: usize| -> [f64; 3] { + let t = through[k]; + let p_in = place(uniform[i], uniform[j], sign).map(|c| inner * c); + let angle = |s: f64| (s * std::f64::consts::FRAC_PI_4).tan(); + let d = place(angle(clustered[i]), angle(clustered[j]), sign); + let len = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt(); + let p_out = d.map(|c| radius * c / len); + [0, 1, 2].map(|c| (1.0 - t) * p_in[c] + t * p_out[c]) + }; + for k in 0..layers { + for j in 0..n { + for i in 0..n { + hexes.push([ + point(i, j, k), + point(i + 1, j, k), + point(i + 1, j + 1, k), + point(i, j + 1, k), + point(i, j, k + 1), + point(i + 1, j, k + 1), + point(i + 1, j + 1, k + 1), + point(i, j + 1, k + 1), + ]); + } + } + } + } + } + let mut elements = Vec::new(); + for mut h in hexes.into_iter().filter(|h| kept(h)) { + // Positive volume: the bottom face counter-clockwise seen from the top. + let sub = |a: [f64; 3], b: [f64; 3]| [a[0] - b[0], a[1] - b[1], a[2] - b[2]]; + let (a, b, c) = (sub(h[1], h[0]), sub(h[3], h[0]), sub(h[4], h[0])); + let det = a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0]); + if det < 0.0 { + h = [h[4], h[5], h[6], h[7], h[0], h[1], h[2], h[3]]; + } + let nodes: Vec = h.iter().map(|x| node(m, *x)).collect(); + m.model.set_element(m.next_element, ElementKind::Hex8, nodes).unwrap(); + elements.push(m.next_element); + m.next_element += 1; + } + let pole = [0.0, 0.0, towards]; + let exterior = m.model.exterior_faces(); + let faces: Vec = m + .faces_where(&elements, |c| { + let r = [c[0] - center[0], c[1] - center[1], c[2] - center[2]]; + let len = (r[0] * r[0] + r[1] * r[1] + r[2] * r[2]).sqrt(); + len > 0.9 * radius && (r[0] * pole[0] + r[1] * pole[1] + r[2] * pole[2]) / len > (std::f64::consts::PI / 6.0).cos() + }) + .into_iter() + .filter(|f| exterior.binary_search(f).is_ok()) + .collect(); + (elements, faces) +} + +/// Hertz's impact of two identical elastic spheres (radius R, mass m, E, ν) at relative speed v: +/// with R* = R/2, m* = m/2, E* = E / (2(1 − ν²)), the approach reaches +/// δ = (15 m* v² / (16 E* √R*))^{2/5}, the force F = (4/3) E* √R* δ^{3/2}, and the contact lasts +/// 2.9432 δ / v (K. L. Johnson, *Contact Mechanics*, Cambridge 1985, §11.4). +fn hertz(radius: f64, young: f64, poisson: f64, density: f64, v: f64) -> (f64, f64, f64) { + let mass = 4.0 / 3.0 * std::f64::consts::PI * radius.powi(3) * density; + let (r_eff, m_eff, e_eff) = (radius / 2.0, mass / 2.0, young / (2.0 * (1.0 - poisson * poisson))); + let delta = (15.0 * m_eff * v * v / (16.0 * e_eff * r_eff.sqrt())).powf(0.4); + (delta, 4.0 / 3.0 * e_eff * r_eff.sqrt() * delta.powf(1.5), 2.9432 * delta / v) +} + +/// Two spheres in a head-on collision (a quarter of each, between the symmetry planes x = 0 and +/// y = 0), with a contact pair between their facing caps: the peak force and the contact +/// duration against Hertz. `n` and `layers` set the mesh; `a_over_r` the contact radius the speed +/// is chosen for. +fn sphere_impact(n: usize, layers: usize, lateral: f64, radial: f64, a_over_r: f64) -> (f64, f64, f64, f64, ExplicitResult) { + let (radius, young, poisson, density) = (1.0, 1.0, 0.3, 1.0); + // The speed that gives a contact radius a = √(R* δ). + let delta = a_over_r * a_over_r * radius * radius / (radius / 2.0); + let mass = 4.0 / 3.0 * std::f64::consts::PI * density; + let e_eff = young / (2.0 * (1.0 - poisson * poisson)); + let v = (delta.powf(2.5) * 16.0 * e_eff * (radius / 2.0f64).sqrt() / (15.0 * mass / 2.0)).sqrt(); + let (delta, force, duration) = hertz(radius, young, poisson, density, v); + let mut m = Mesh::new(); + let (upper, upper_cap) = quarter_ball(&mut m, [0.0, 0.0, radius], radius, n, layers, lateral, radial, -1.0); + let (lower, lower_cap) = quarter_ball(&mut m, [0.0, 0.0, -radius], radius, n, layers, lateral, radial, 1.0); + m.material("Mat", young, poisson, density); + m.section("upper", &upper, "Mat"); + m.section("lower", &lower, "Mat"); + m.face_set("upper_cap", upper_cap); + m.face_set("lower_cap", lower_cap); + m.pair("impact", "lower_cap", "upper_cap", Friction::NONE); + let both: Vec = upper.iter().chain(&lower).copied().collect(); + m.node_set("x0", m.nodes_where(&both, |x| x[0].abs() < 1e-9)); + m.node_set("y0", m.nodes_where(&both, |x| x[1].abs() < 1e-9)); + let held = |set: &str, k: usize| { + let mut values = [None; 6]; + values[k] = Some(0.0); + BoundaryCondition::Displacement { target: Target::Set(set.into()), values } + }; + let mut step = ExplicitDynamic::new(1.15 * duration); + step.initial_velocities = vec![ + InitialVelocity { set: "upper".into(), velocity: [0.0, 0.0, -v / 2.0] }, + InitialVelocity { set: "lower".into(), velocity: [0.0, 0.0, v / 2.0] }, + ]; + step.bulk_viscosity = BulkViscosity { linear: 0.0, quadratic: 1.5 }; + step.history_interval = Some(duration / 2000.0); + step.output_interval = Some(1.15 * duration); + m.step(step, &[("x0", held("x0", 0)), ("y0", held("y0", 1))], &[]); + let r = m.solve(); + let s = &r.contacts[0]; + // A quarter of each sphere carries a quarter of the force. + let peak = 4.0 * s.peak_force; + let lasted = s.last_contact.unwrap_or(0.0) - s.first_contact.unwrap_or(0.0); + eprintln!( + "Hertz ({} elements, a/R = {a_over_r}): v = {v:.4e}, δ = {delta:.4e}; peak force {peak:.5e} vs {force:.5e} ({:+.2}%), duration {lasted:.4} vs {duration:.4} ({:+.2}%); max penetration {:.2e}; KE after/before {:.4}; balance {:.2e}; {} increments", + both.len(), + 100.0 * (peak / force - 1.0), + 100.0 * (lasted / duration - 1.0), + s.max_penetration, + r.energy.kinetic / r.initial_kinetic, + r.energy.relative_error(r.initial_kinetic), + r.increments + ); + (peak, force, lasted, duration, r) +} + +/// A coarse sphere impact in the main gate: Hertz's force and duration to a few percent. +#[test] +fn spheres_collide_as_hertz_says_coarse() { + let (peak, force, lasted, duration, _) = sphere_impact(12, 6, 1.25, 1.25, 0.15); + assert!((peak - force).abs() < 0.1 * force, "{peak} vs {force}"); + assert!((lasted - duration).abs() < 0.1 * duration, "{lasted} vs {duration}"); +} + +/// The fine sphere impact: contact radius a tenth of the radius, about six elements across it. +#[test] +#[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"] +fn spheres_collide_as_hertz_says() { + let (peak, force, lasted, duration, r) = sphere_impact(24, 12, 1.23, 1.185, 0.1); + assert!((peak - force).abs() < 0.03 * force, "{peak} vs {force}"); + assert!((lasted - duration).abs() < 0.03 * duration, "{lasted} vs {duration}"); + assert!(r.energy.kinetic > 0.95 * r.initial_kinetic); +} diff --git a/crates/explicit/tests/large_deformation.rs b/crates/explicit/tests/large_deformation.rs index 5901806..c195b16 100644 --- a/crates/explicit/tests/large_deformation.rs +++ b/crates/explicit/tests/large_deformation.rs @@ -168,7 +168,7 @@ impl Driver { let resolved = model.resolve().unwrap(); let step = model.step("Dynamic").unwrap(); let Procedure::ExplicitDynamic(settings) = &step.procedure else { panic!() }; - let sys = System::new(model, &resolved, step, settings, 1000).unwrap(); + let sys = System::new(model, &resolved, step, settings, &Limits { max_dofs: 1000, ..Limits::default() }).unwrap(); let (state, ws) = (State::new(&sys), Workspace::new(&sys)); let u = vec![0.0; sys.dofs()]; Driver { sys, state, ws, u, time: 0.0 } diff --git a/crates/fem/src/contact.rs b/crates/fem/src/contact.rs index 714cafa..971868c 100644 --- a/crates/fem/src/contact.rs +++ b/crates/fem/src/contact.rs @@ -311,6 +311,32 @@ impl Grid { } } + /// The boxes that may overlap box `b`: those listed in the cells it overlaps, in box order + /// without repeats (into `out`, cleared first). + pub fn candidates_in(&self, b: &Aabb, out: &mut Vec) { + out.clear(); + if self.items.is_empty() || b.iter().flatten().any(|c| !c.is_finite()) { + return; + } + let top = [0, 1, 2].map(|k| self.lowest[k] + self.cell * self.dims[k] as f64); + if (0..3).any(|k| b[1][k] < self.lowest[k] || b[0][k] > top[k]) { + return; + } + let (lo, hi) = (self.coordinates(b[0]), self.coordinates(b[1])); + for k in lo[2]..=hi[2] { + for j in lo[1]..=hi[1] { + for i in lo[0]..=hi[0] { + let c = self.index([i, j, k]); + if let (Some(a), Some(e)) = (self.start.get(c), self.start.get(c + 1)) { + out.extend_from_slice(self.items.get(*a..*e).unwrap_or_default()); + } + } + } + } + out.sort_unstable(); + out.dedup(); + } + /// The number of entries (box × cell overlaps). pub fn entries(&self) -> usize { self.items.len() @@ -394,5 +420,12 @@ mod tests { assert!(g.entries() <= MAX_GRID_ITEMS); } assert_eq!(Grid::new(&[], 10).candidates([0.0; 3]), &[] as &[u32]); + // Box queries: every box that overlaps, once, in order. + let mut found = Vec::new(); + grid.candidates_in(&[[10.2, 0.2, 0.2], [12.1, 0.4, 0.4]], &mut found); + let expected: Vec = (0..50).filter(|i| boxes[*i as usize][0][0] <= 12.1 && boxes[*i as usize][1][0] >= 10.2).collect(); + assert!(expected.iter().all(|e| found.contains(e)) && found.windows(2).all(|w| w[0] < w[1]), "{found:?} {expected:?}"); + grid.candidates_in(&[[-9.0, 0.0, 0.0], [-8.0, 1.0, 1.0]], &mut found); + assert!(found.is_empty()); } } diff --git a/crates/testkit/src/oracle.rs b/crates/testkit/src/oracle.rs index f1cc447..c51e773 100644 --- a/crates/testkit/src/oracle.rs +++ b/crates/testkit/src/oracle.rs @@ -263,7 +263,8 @@ pub fn run_case(dir: &Path, expected: &Expected) -> Result, m: &Value) -> Result { let pointer = m.get("pointer").and_then(Value::as_str).ok_or("no `pointer`")?; let fresh; @@ -278,7 +279,15 @@ fn measure(session: &mut Session, outputs: &BTreeMap, m: &Value) }; let v = out.pointer(pointer).ok_or_else(|| format!("nothing at `{pointer}`"))?; let x = v.as_f64().ok_or_else(|| format!("`{pointer}` is {v}, not a number"))?; - let x = x * m.get("scale").and_then(Value::as_f64).unwrap_or(1.0); + let mut x = x * m.get("scale").and_then(Value::as_f64).unwrap_or(1.0); + // A ratio: this value over another measure's (both from the same run). + if let Some(over) = m.get("over") { + let y = measure(session, outputs, over)?; + if y == 0.0 { + return Err(format!("`{pointer}` is measured over a zero")); + } + x /= y; + } if x.is_finite() { Ok(x) } else { Err(format!("`{pointer}` is not finite")) } } @@ -730,6 +739,18 @@ mod tests { assert_eq!(fixture_band(1e-16, 0.0, a, FIXTURE_RELATIVE), 1e-9, "a zero reference keeps its tolerance as the band"); } + #[test] + fn measures_read_pointers_scale_and_divide() { + let mut session = Session::default(); + let outputs: BTreeMap = [("s".to_string(), json!({"a": 6.0, "b": -2.0, "z": 0.0}))].into_iter().collect(); + let m = |v: Value| measure(&mut Session::default(), &outputs, &v); + assert_eq!(m(json!({"from": "s", "pointer": "/a", "scale": 0.5})), Ok(3.0)); + assert_eq!(m(json!({"from": "s", "pointer": "/a", "over": {"from": "s", "pointer": "/b", "scale": -1.0}})), Ok(3.0)); + assert!(m(json!({"from": "s", "pointer": "/a", "over": {"from": "s", "pointer": "/z"}})).is_err()); + assert!(m(json!({"from": "s", "pointer": "/a", "over": {"from": "s", "pointer": "/nope"}})).is_err()); + assert!(measure(&mut session, &outputs, &json!({"from": "t", "pointer": "/a"})).is_err()); + } + #[test] fn expected_files_are_checked() { let good = json!({"title": "t", "kind": "theory", "source": "s", "quantities": [ diff --git a/docs/architecture.md b/docs/architecture.md index 7c116bf..ea47dda 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -255,7 +255,7 @@ number their degrees of freedom with `fieldcraft_fem::DofLayout` (docs/fem.md). increment. Other element kinds and formulations, couplings, reference points, rotary inertia, rotational springs and moments return "not supported yet". Masses are row-sum lumped. A zero-valued boundary condition holds a degree of freedom; a displacement value is reached at constant velocity over the step. Loads act in full from the - start. The step lists initial velocities and fixed frictionless rigid planes by set. + start. The step lists initial velocities and fixed rigid planes (with friction) by set. 2. **Time step.** Each element's highest frequency is bounded in closed form (exact for a cube); the increment is the step's scale (0.9 by default) times the smallest element increment with Wilkins' bulk-viscosity correction, recomputed every increment. The ten controlling elements @@ -280,6 +280,15 @@ number their degrees of freedom with `fieldcraft_fem::DofLayout` (docs/fem.md). `fem::plasticity`. Elements whose damage reaches 1 or whose plastic strain reaches the material's erosion strain are eroded (reported with their energy; their nodes keep mass and momentum). The energy account adds the plastic work and the eroded energy. +7. **Contact (X3).** The model's contacts (`interaction.contact.*`) become interfaces before the + elements are built (nodes that start inside a surface are moved onto it or refused): + segment-to-segment penalty contact both ways (projection, polygon clipping and quadrature over + each pair of facing faces), a bucket-grid search every increment in parallel by faces, a + penalty per area from the masses and the stable increment that never shortens it, and Coulomb + friction by a return map at the secondary nodes; rigid walls take Coulomb impulses. The contact + forces join the equations of motion beside the external ones, their work and the friction + dissipation join the energy account, and each contact's peak force, first and last contact, + deepest penetration and friction are reported (docs/fem.md). Verification (`crates/explicit/tests/verification.rs` and `large_deformation.rs`, analytic only): stress waves in hex8 and tet4 bars (front at c·t, σ = ρ c v, work input), a bar striking a @@ -288,7 +297,10 @@ bit-for-bit repeatability; Dienes' simple shear for both rates, a stressed eleme by 90°, the J2 uniaxial curve and unloading at large strain, damage failure at the fracture strain, and the Taylor bar against Taylor's 1948 estimate; `crates/explicit/tests/one_d.rs`: spring–mass period, damped decay rate and period, a dashpot alone, two masses on a line spring, -waves in rod2 and rod3 bars; in the engine (`solve.rs`), a +waves in rod2 and rod3 bars; `crates/explicit/tests/contact.rs`: colliding bars, sliding +blocks under constant and decaying friction on a rigid wall and on a plate, the contact patch +test, Hertz's sphere impact, a free body untouched by general contact, initial penetrations and +thread-count independence; in the engine (`solve.rs`), a cantilever's free-vibration period against `solve.frequency`. ## Lattice Boltzmann solver diff --git a/docs/fem.md b/docs/fem.md index 809dfdf..655438a 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -347,3 +347,59 @@ form (exact), and a tube inflated by a follower pressure while stretched axially block at a thermal stretch of 1.2 (exact, full and hybrid), a bar held at its ends heated past yield (exact), the linear limit against the linear and coupled solvers, and tangents with a temperature field. + +crushing stops at the minimum increment, naming it; the caps; what is refused). + +## Contact and friction (explicit dynamics) + +**Model** (`fieldcraft_model::contact`): a contact names its surfaces as sets of element faces — +a **pair** (`main`, `secondary`) or **general** contact on one surface (self contact) or on every +exterior face of the solids — with Coulomb **friction** μ = μ_k + (μ_s − μ_k) e^{−d γ̇} (static, +kinetic, decay per slip rate), a **penalty scale** s ≤ 1 and what happens to nodes that start +inside the other surface (**adjust**: moved onto it strain-free and reported; **refuse**: the +analysis stops naming the node, the face and the depth). Planar rigid walls take the same friction. +Commands `interaction.contact.pair|general|edit|delete|find`; Abaqus writes `*SURFACE +INTERACTION` with `*FRICTION` (exponential decay), `*CONTACT PAIR` and general `*CONTACT`, LS-DYNA +`*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE` and `_SINGLE_SURFACE` (`fs`, `fd`, `dc`, `sfsa`, +`sfsb`); Nastran and the implicit solver refuse contact by name (implicit contact is M7). + +**Geometry** (`fieldcraft_fem::contact`, shared with the implicit solver to come): the closest +point of a triangle (exact region tests) or of a bilinear quadrilateral (projected Newton with the +exact Hessian) to a point, with the signed gap along the face's normal and the lateral offset; +a bucket grid of face boxes whose cell is twice the mean box, grown until the grid and its entries +fit their caps, queried by point or by box with the candidates in face order. + +**Kinematics** (`fieldcraft_explicit::contact`): segment to segment. Each secondary face and +each main face facing it (normals within 84° of opposite) are projected onto the secondary face's +plane, split into triangles and clipped against each other; on each overlap triangle a +three-point rule takes the gap g = (x_s − x_m)·n_m between the faces' closest points. Where +g < 0 (no deeper than half the main face's element), a penalty pressure p = ε(−g) acts on the +secondary face's nodes with weights N_i^s along n_m and on the main face's nodes with −N_j^m. For +a uniform pressure across flat faces both meshes receive exactly the consistent nodal forces, +so the contact patch test passes on non-matching meshes. Every contact runs both ways (a pair +with each surface as the secondary, single-surface contact with every face against every other +face that shares no node with it) at half the penalty each. + +**Penalty**: ε = s γ min_i(m_i/A_i)/Δt² over the contact's nodes (m_i lumped mass, A_i the node's +area facing one way, Δt the stable increment, γ = 0.3). By Gershgorin's and Weyl's inequalities +the contact adds at most 2 s γ/Δt² to the square of any frequency, which the default increment +scale 0.9 leaves room for: contact never shortens the increment. + +**Friction**: at the secondary nodes, by a regularised return map against the main face that +pressed each node most: the elastic predictor t* = P t_n − k_t Δu_T (k_t = ε A_i/2), stick when +|t*| ≤ μ f_n, otherwise slip back to the cone, dissipating μ f_n (|t*| − μ f_n)/k_t. On rigid +walls, a node's tangential velocity changes by at most μ times its normal velocity change. + +**Energy**: the contact forces' work (penalty energy plus friction dissipated) is integrated with +the trapezoidal rule and enters the balance; friction dissipation (contact and walls) is +reported apart. **Search** every increment, in parallel by faces, forces added in face order: the +results are the same on any number of threads. + +**Verification** (`crates/explicit/tests/contact.rs`, `tests/oracle/explicit-*`): two bars +colliding (contact for 2L/c, interface force ρcAv₀/2, velocities exchanged, kinetic energy kept to +0.01%, the contact energy back to 10⁻⁴ of it), a block sliding to a stop under constant and +decaying friction on a rigid wall and on a deformable plate (distance and time to 0.1–0.5%), the +contact patch test on non-matching meshes (σ_zz uniform to 0.1%), two spheres colliding against +Hertz (peak force +2.3%, duration −0.4% on a graded quarter-sphere mesh), a free body untouched +by general contact on all its faces, initial penetrations adjusted or refused, the face cap, and +bit-for-bit results on one and four threads. diff --git a/docs/oracle.md b/docs/oracle.md index e66ad67..def31ef 100644 --- a/docs/oracle.md +++ b/docs/oracle.md @@ -2,7 +2,7 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the policy is [verification.md](verification.md). Every case is our own model, built with FieldCraft commands from the problem's published description, and compared with the reference value its `expected.json` cites, within the stated tolerance and for the stated reason. Accepted results are frozen in each case's `fixture.json`: a later change that moves one beyond its band (10⁻⁹ of its size) fails CI even while the reference still passes. No reference value was produced with a commercial product. -**176 of 182 quantities pass** (6 known failure(s)) in 53 active cases: 44 theory, 9 published benchmark, 0 experiment; 9 pending cases wait for capabilities still being built. +**193 of 199 quantities pass** (6 known failure(s)) in 57 active cases: 48 theory, 9 published benchmark, 0 experiment; 10 pending cases wait for capabilities still being built. ## Theory: closed-form solutions, patch tests and mesh invariants @@ -28,6 +28,23 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the | [cdt-stepped-block](../tests/oracle/cdt-stepped-block/) | volume of the tetrahedra | 47000 | 47000 | < 1e-7% | 1e-12 rel. | pass | matches | | [cdt-stepped-block](../tests/oracle/cdt-stepped-block/) | elements that fail the certified Jacobian check | 0 | 0 | ≈ 0 | exact | pass | matches | | [cdt-stepped-block](../tests/oracle/cdt-stepped-block/) | elements the Jacobian check could not decide | 0 | 0 | ≈ 0 | exact | pass | matches | +| [explicit-bars-collision](../tests/oracle/explicit-bars-collision/) | time of the last contact (the bars touch from t = 0) against 2L/c = 0.2 | 0.2 | 0.20242 | 1.21% | 2% | pass | matches | +| [explicit-bars-collision](../tests/oracle/explicit-bars-collision/) | the largest force across the interface against ρcAv₀/2 = 0.05 | 0.05 | 0.05147 | 2.94% | 4% | pass | matches | +| [explicit-bars-collision](../tests/oracle/explicit-bars-collision/) | velocity of the struck bar's far end after the impact (0: the striker stops) | 0 | -0.0033987 | -3.4e-3 | ± 0.005 | pass | matches | +| [explicit-bars-collision](../tests/oracle/explicit-bars-collision/) | velocity of the struck bar's far end after the impact (v₀ = 0.1) | 0.1 | 0.1034 | 3.4% | 5% | pass | matches | +| [explicit-bars-collision](../tests/oracle/explicit-bars-collision/) | kinetic energy after the impact over ½mv₀² (frictionless elastic impact) | 1 | 0.99995 | -0.00468% | 1% | pass | matches | +| [explicit-bars-collision](../tests/oracle/explicit-bars-collision/) | energy left in the contact after the separation, over ½mv₀² (its drift) | 0 | -7.7435e-5 | -7.7e-5 | ± 0.001 | pass | matches | +| [explicit-contact-patch](../tests/oracle/explicit-contact-patch/) | σ_zz of the lower block's corner element under the interface (element 10) over that of the upper block's corner element over it (element 19) | 1 | 1.0001 | 0.0118% | 0.5% | pass | matches | +| [explicit-contact-patch](../tests/oracle/explicit-contact-patch/) | σ_zz of the lower block's centre element under the interface (element 14) over that of element 19 | 1 | 1.0003 | 0.0287% | 0.5% | pass | matches | +| [explicit-contact-patch](../tests/oracle/explicit-contact-patch/) | σ_zz of the lower block's bottom corner element (element 1) over that of element 19 | 1 | 1.0001 | 0.0134% | 0.5% | pass | matches | +| [explicit-contact-patch](../tests/oracle/explicit-contact-patch/) | σ_zz of the upper block's top far corner element (element 26) over that of element 19 | 1 | 0.99944 | -0.0564% | 0.5% | pass | matches | +| [explicit-contact-patch](../tests/oracle/explicit-contact-patch/) | the force across the interface over −σ_zz A (A = 1) | 1 | 1.0026 | 0.259% | 1% | pass | matches | +| [explicit-contact-patch](../tests/oracle/explicit-contact-patch/) | σ_xx of element 14 over its σ_zz (uniaxial stress: 0) | 0 | -3.6052e-4 | -3.6e-4 | ± 0.005 | pass | matches | +| [explicit-sliding-block-contact](../tests/oracle/explicit-sliding-block-contact/) | distance travelled (top centre node) against ∫ v dv / (g μ(v)) = 2.01524 | 2.0152 | 2.0118 | -0.172% | 1% | pass | matches | +| [explicit-sliding-block-contact](../tests/oracle/explicit-sliding-block-contact/) | velocity of the top centre node after the stop | 0 | ≈ 0 | ≈ 0 | ± 0.03 | pass | matches | +| [explicit-sliding-block-wall](../tests/oracle/explicit-sliding-block-wall/) | distance travelled (top centre node) against v₀²/(2μg) = 1.5 | 1.5 | 1.4999 | -0.0053% | 1% | pass | matches | +| [explicit-sliding-block-wall](../tests/oracle/explicit-sliding-block-wall/) | velocity of the top centre node after the stop | 0 | ≈ 0 | ≈ 0 | ± 0.03 | pass | matches | +| [explicit-sliding-block-wall](../tests/oracle/explicit-sliding-block-wall/) | friction dissipation over ½mv₀² = 2.25 | 1 | 1.0012 | 0.121% | 1% | pass | matches | | [filleted-bracket-fillet-stress](../tests/oracle/filleted-bracket-fillet-stress/) | largest nodal von Mises stress on the fillet face (MPa) | 213 | mesh-dependent | within 3% | 3% | pass | — | | [filleted-bracket-fillet-stress](../tests/oracle/filleted-bracket-fillet-stress/) | total reaction x (N) | -2000 | -2000 | < 1e-7% | 1e-9 rel. | pass | matches | | [filleted-bracket-fillet-stress-coarse](../tests/oracle/filleted-bracket-fillet-stress-coarse/) | largest nodal von Mises stress on the fillet face (MPa) | 213 | mesh-dependent | within 3% | 3% | pass | — | @@ -206,6 +223,10 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the - **[cantilever-tet10](../tests/oracle/cantilever-tet10/)** (theory): Cantilever under a tip load, Timoshenko deflection (tet10). Source: Timoshenko and Gere, Mechanics of Materials (beam with shear deformation: δ = PL³/3EI + PL/κGA, κ = 5/6 for a rectangle); Cowper, J. Appl. Mech. 33 (1966) 335–340 for the shear coefficient. Tolerances — tip_uz: beam theory vs 3D elasticity (a few tenths of a percent at L/h = 10); quadratic tetrahedra are converged on this mesh; reaction: equilibrium holds to rounding - **[cdt-cylinder](../tests/oracle/cdt-cylinder/)** (theory): Constrained Delaunay mesh of a cylinder: volume within the chordal bound, valid elements. Source: mesh invariants: an inscribed polygon of sagitta s loses at most 2s/R of a circle's area (elementary geometry); constrained Delaunay tetrahedralization after Diazzi et al., ACM TOG 42(6) 2023. Tolerances — volume: the tessellation inscribes the cylinder within the chordal tolerance s = 0.05: the polygon section loses between 0 and 2s/R = 0.5% of the area; the reference is the middle of that band; invalid: every element of a valid mesh has positive volume - **[cdt-stepped-block](../tests/oracle/cdt-stepped-block/)** (theory): Constrained Delaunay mesh of a stepped block: exact volume, valid elements. Source: mesh invariants: a tetrahedralization of a polyhedron fills its volume exactly (Si, TetGen paper ACM TOMS 41(2) 2015 for the definition of a CDT; Diazzi et al., ACM TOG 42(6) 2023); the volume is 60·40·10 + 40·40·10 + 20·40·10 − 10·10·10. Tolerances — volume: a planar boundary is tessellated exactly and the constrained Delaunay mesh fills it exactly: only rounding; invalid: every element of a valid mesh has positive volume; undetermined: a linear tetrahedron's sign is exact +- **[explicit-bars-collision](../tests/oracle/explicit-bars-collision/)** (theory): Two identical elastic bars colliding (explicit contact). Source: One-dimensional elastic impact of equal bars: contact for 2L/c, interface force ρcAv₀/2, exchange of velocities, no loss of energy; Graff, Wave Motion in Elastic Solids, Dover 1991, §2.1–2.3; Johnson, Impact Strength of Materials, Arnold 1972, §1.6. Tolerances — duration: the penalty springs at the interface delay the release by a few element transit times (about 1% of 2L/c with 200 elements per bar); force: the plateau is ρcAv₀/2 to 0.01%; the peak is the plateau plus the ringing behind the discrete wave front, a few percent; v_a: the striker stops on average (its kinetic energy is under 0.01% of the start's, see ke); its far end node keeps the ringing of the discrete release front, a few percent of v₀; v_b: the struck bar leaves at v₀ on average; its far end node keeps the ringing of the discrete wave fronts, a few percent; ke: an elastic impact keeps the kinetic energy; the quadratic bulk viscosity at the compression fronts and the vibration left in the bars take a fraction of a percent; contact_energy: the penalty springs give back what they took; the trapezoidal work of a spring that engages between increments drifts by far less than 0.1% +- **[explicit-contact-patch](../tests/oracle/explicit-contact-patch/)** (theory): Contact patch test: two blocks with non-matching meshes in uniaxial compression (explicit contact). Source: The contact patch test: a contact discretisation must transmit a uniform pressure across a non-matching interface exactly; Taylor & Papadopoulos, On a patch test for contact problems in two dimensions, Nonlinear Computational Mechanics, Springer 1991, 690–702; Crisfield, Re-visiting the contact patch test, IJNME 48:435–449, 2000. Uniaxial stress between frictionless symmetry planes is uniform in both blocks. Tolerances — lower_corner: a uniform stress is an equilibrium of the discrete system (the integrated contact forces are the consistent nodal forces of a uniform pressure on both meshes); what is left is the ringing of the loading, a hundred wave transits slow; lower_centre: as for the corner element; lower_base: as for the corner element; upper_top: as for the corner element; force: the integrated contact pressure carries the whole load; the largest force of the run is the final one in a slow squeeze, up to the ringing; lateral: frictionless symmetry planes and a frictionless interface leave the blocks free to expand sideways; the non-matching interface must not hold them +- **[explicit-sliding-block-contact](../tests/oracle/explicit-sliding-block-contact/)** (theory): Block sliding to a stop on a deformable plate, friction decaying with the slip rate (explicit contact). Source: A rigid block launched at v₀ with friction μ(v) = μk + (μs − μk) e^(−d v) decelerates at g μ(v); it travels ∫₀^v₀ v dv / (g μ(v)) (Simpson's rule, 200 000 intervals) and stops after [v₀ + ln((μk + (μs − μk) e^(−d v₀))/μs)/d]/(g μk) = 1.2001; Meriam & Kraige, Engineering Mechanics: Dynamics, §3/5; the exponential decay law as in Oden & Martins, CMAME 52:527–634, 1985. Tolerances — distance: the block bounces on the penalty springs and the plate under it shears elastically, so the normal force and the slip rate it feels vary by a little around their rigid-body values; velocity: the block sticks; what is left is its elastic vibration, under 1% of v₀ +- **[explicit-sliding-block-wall](../tests/oracle/explicit-sliding-block-wall/)** (theory): Block sliding to a stop on a rigid wall with Coulomb friction (explicit). Source: A rigid block launched at v₀ on a horizontal plane with Coulomb friction μ decelerates at μg and stops after v₀/(μg) at v₀²/(2μg), its kinetic energy all dissipated by friction; Meriam & Kraige, Engineering Mechanics: Dynamics, §3/5 (work and energy). Tolerances — distance: the wall's friction impulses are exact Coulomb impulses; the block's elastic vibration from the sudden gravity and the friction moment move the node by a fraction of a percent; velocity: the block sticks; what is left is its elastic vibration, under 1% of v₀; dissipated: the kinetic energy goes into friction; the vibration and the work of gravity as the block settles add or take a fraction of a percent - **[filleted-bracket-fillet-stress](../tests/oracle/filleted-bracket-fillet-stress/)** (theory): Filleted L-bracket under a side load: the peak von Mises stress on the R6 fillet, fine tet10 mesh. Source: mesh convergence of the finite element solution (P. J. Roache, Verification and Validation in Computational Science and Engineering, Hermosa, 1998, ch. 5): the bracket has no closed-form stress; the reference is the mesh-independent peak of a four-level study of quadratic tetrahedral meshes refined by √2 in size and in curvature resolution (README.md); equilibrium. Tolerances — fillet_peak: the four-level study gives 213.3, 210.6, 214.3, 213.3 MPa (2k to 13.5k elements): the bending stress is converged across the fillet, and the peak nodal value moves only with which of the fillet's tessellation lines (fixed angles, through every level) carries the largest node, by under 2%; reaction: equilibrium holds to rounding - **[filleted-bracket-fillet-stress-coarse](../tests/oracle/filleted-bracket-fillet-stress-coarse/)** (theory): Filleted L-bracket under a side load: the peak von Mises stress on the R6 fillet, coarse tet10 mesh. Source: mesh convergence of the finite element solution (P. J. Roache, Verification and Validation in Computational Science and Engineering, Hermosa, 1998, ch. 5): the bracket has no closed-form stress; the reference is the mesh-independent peak of a four-level study of quadratic tetrahedral meshes refined by √2 in size and in curvature resolution (README.md); equilibrium. Tolerances — fillet_peak: the four-level study gives 213.3, 210.6, 214.3, 213.3 MPa (2k to 13.5k elements): the bending stress is converged across the fillet, and the peak nodal value moves only with which of the fillet's tessellation lines (fixed angles, through every level) carries the largest node, by under 2%; reaction: equilibrium holds to rounding - **[free-beam-frequencies](../tests/oracle/free-beam-frequencies/)** (theory): Natural frequencies of a free beam: rigid-body modes and the first bending mode. Source: Euler–Bernoulli free–free beam, β₁L = 4.7300; Timoshenko, Young and Weaver, Vibration Problems in Engineering, 4th ed., Wiley 1974, §5.3. Tolerances — rigid_6: a free body has six zero frequencies; numerical rigid-body frequencies below 10⁻³ of the first elastic one (eigenvalues below 10⁻⁶ of it) are rounding of the shifted solve; f7: Euler–Bernoulli neglects shear and rotary inertia, about 0.8% at L/h = 20 for the free–free mode; mass: the mass matrix integrates the density exactly @@ -259,6 +280,7 @@ Reference values already cited; each case gets its recipe when the capability it | Case | Kind | Problem | Needs | Reference values | Source | |---|---|---|---|---|---| | [cylinder-channel-2d1](../tests/oracle/cylinder-channel-2d1/) | benchmark | Flow around a cylinder, steady 2D-1 (Re = 20) | lattice Boltzmann CFD (M12) | drag coefficient (published interval 5.57–5.59) = 5.58; lift coefficient (interval 0.0104–0.0110) = 0.0107 | Schäfer and Turek, 'Benchmark computations of laminar flow around a cylinder', Notes on Numerical Fluid Mechanics 52 (1996) 547–566, Table 4 (reference intervals) | +| [explicit-hertz-spheres](../tests/oracle/explicit-hertz-spheres/) | theory | Two elastic spheres in a head-on collision: Hertz's force and contact duration (explicit contact) | a graded sphere mesh from commands (a cubed-sphere hex mesher, or tet refinement towards a point); until then the explicit crate's tests `spheres_collide_as_hertz_says` (a graded cubed-sphere quarter of each sphere, 27 648 hex8: peak force +2.3%, duration −0.4%) and its coarse twin in the main gate check this case on a mesh generated in code | the largest contact force, two steel spheres R = 10 mm (E = 210 GPa, ν = 0.3, ρ = 7.85 t/m³) meeting at 1 m/s (N) = 1005.4; the contact duration (s) = 6.0163e-5 | Hertz's quasi-static impact: with R* = R/2, m* = m/2, E* = E/(2(1 − ν²)), the approach reaches δ = (15 m* v²/(16 E* √R*))^(2/5), the force (4/3) E* √R* δ^(3/2), and the contact lasts 2.9432 δ/v; Hertz, J. reine angew. Math. 92:156–171, 1882; K. L. Johnson, Contact Mechanics, Cambridge 1985, §11.4 | | [lid-driven-cavity-re100](../tests/oracle/lid-driven-cavity-re100/) | benchmark | Lid-driven cavity, Re = 100 | lattice Boltzmann CFD (M12) | u/U at y = 0.4531 on the vertical centreline = -0.2109 | Ghia, Ghia and Shin, Journal of Computational Physics 48 (1982) 387–411, Table I (u along the vertical centreline) | | [mh-twisted-beam](../tests/oracle/mh-twisted-beam/) | benchmark | MacNeal–Harder twisted beam | shell elements (M8) or mapped blocks with helical edges | tip displacement along the in-plane load = 0.005424; tip displacement along the out-of-plane load = 0.001754 | MacNeal and Harder, 'A proposed standard set of problems to test finite element accuracy', Finite Elements in Analysis and Design 1 (1985) 3–20, Fig. 6 (p. 8), Table 3 (p. 10) | | [nafems-le11](../tests/oracle/nafems-le11/) | benchmark | NAFEMS LE11 solid cylinder/taper/sphere under a temperature field | thermal loads (M8 heat transfer) | σ_zz at A = -105 | NAFEMS, The Standard NAFEMS Benchmarks, TNSB Rev. 3, October 1990, test LE11 (target σ_zz = −105 MPa at A) | diff --git a/docs/parity.md b/docs/parity.md index c56fff9..d7b4719 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -55,7 +55,7 @@ products are run, no documentation text is copied). Status: **done**, **started* | Model | Steps: linear static, frequency; output presets | M3 | | done | | Model | Steps: nonlinear static, explicit dynamics, heat transfer; amplitudes | M7–M8, M15 | | started (nonlinear static with or without nonlinear geometry, increments, tolerances and a load path (cyclic loading), written to Abaqus `*STEP, NLGEOM` and LS-DYNA implicit (the load path as LS-DYNA's curve; refused by Abaqus for now), refused by Nastran; explicit dynamics; heat transfer; amplitudes per load planned) | | Interaction | Tie between non-matching meshes | M7 | | planned | -| Interaction | Contact (general and pairs), friction, automatic pair detection | M7 | | planned | +| Interaction | Contact (general and pairs), friction, automatic pair detection | M7, M10 | | started (explicit: pairs, general and self contact, segment to segment with a penalty that keeps the time increment, Coulomb friction with exponential decay, friction on rigid walls, `interaction.contact.find`, Abaqus and LS-DYNA contact keywords; implicit contact in M7) | | Interaction | Springs and dashpots (ground or two-point; DOF, vector or line), point masses with rotary inertia | M4 | | done (implicit; explicit translational; all three decks with what each cannot take reported) | | Interaction | Bolt pretension, connectors | M7 | | planned | | Decks | Abaqus input writer (`deck.write`, flat deck; reports what Abaqus cannot express) | M3 | | done | @@ -64,7 +64,7 @@ products are run, no documentation text is copied). Status: **done**, **started* | 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; 2D plane stress, plane strain and axisymmetric tri3, tri6, quad4 (full, incompatible modes), quad8 (full, reduced); quad4/tri3 MITC shells and Timoshenko beams with profiles, mixed with solids; reference points, kinematic and distributing couplings, moments; centrifugal loads) | | Solve | Own nonlinear implicit solver: Newton–Raphson with line search and automatic increments, total Lagrangian solids and 2D elements, mixed u–p at finite strain, follower pressure | M15 | yes | started (solids, plane strain, axisymmetric; St Venant–Kirchhoff, neo-Hookean and J2 plasticity (Simo's exponential return map at finite strain, radial return at small strain, consistent tangents, the return shared with explicit dynamics); verified by the elastica, tube inflation, patch and objectivity tests, uniaxial tension to large strain, Hill's thick cylinder, a cantilever's collapse and the Bauschinger effect; Mooney–Rivlin and Yeoh; thermal strain at finite and small strain; a symmetric indefinite LDLᵀ factorization with inertia, Newton through indefinite tangents (hybrid elements at ν = 0.4999 under prescribed displacement); shells, contact, buckling, arc-length next) | -| 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 and X2: large deformation, objective rates, J2 and Johnson–Cook plasticity, damage, erosion; contact planned) | +| 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 and X2: large deformation, objective rates, J2 and Johnson–Cook plasticity, damage, erosion; X3: deformable contact and friction, verified by bar collision, sliding blocks, the contact patch test and Hertz) | | 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 | | Results | Readers for public result formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`) | M9 | | planned | diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index 10e6f1e..5ff2ecb 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -1246,3 +1246,43 @@ update are written from (cited again in the module docs): 34(1):82–109, 1992. - Iterative refinement: N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM 2002, §12.2. + +## Part Q: Contact and friction (explicit), as implemented + +The papers and books `fieldcraft_fem::contact` and `fieldcraft_explicit::contact` are written +from (cited again in the module docs); docs/fem.md describes the formulation. + +- Closest point on a triangle by its Voronoi regions: C. Ericson, *Real-Time Collision + Detection*, Morgan Kaufmann 2005, §5.1.5; on a bilinear quadrilateral by Newton's method on the + stationarity conditions, projected on the parameter square: the local search of node-to-segment + contact, J. O. Hallquist, G. L. Goudreau, D. J. Benson, CMAME 51:107–137, 1985 (A5.1), and P. + Wriggers, *Computational Contact Mechanics*, 2nd ed., Springer 2006, §4.1. +- Global search by bucket sorting into a uniform grid: D. J. Benson, J. O. Hallquist, CMAME + 78:141–163, 1990 (A5.2). We chose a grid over the BVH of A5.4: one pass per increment, memory + bounded by caps whatever the mesh, candidates in face order (determinism); the cell grows when a + few huge faces would blow the entry cap. +- Segment-to-segment integration (projection on the secondary face's plane, clipping, + sub-triangle quadrature): M. A. Puso, T. A. Laursen, *A mortar segment-to-segment contact + method for large deformation solid mechanics*, CMAME 193:601–629, 2004; with the penalty applied + at the integration points, K. A. Fischer, P. Wriggers, Comput. Mech. 36:226–244, 2005; polygon + clipping, I. E. Sutherland, G. W. Hodgman, Commun. ACM 17(1):32–42, 1974. Chosen over + node-to-segment because the latter fails the contact patch test on non-matching meshes (R. L. + Taylor, P. Papadopoulos, in *Nonlinear Computational Mechanics*, Springer 1991, 690–702; M. A. + Crisfield, IJNME 48:435–449, 2000), which our first node-to-segment version showed too (12% + spread of σ_zz in the patch test, 0.1% now). +- Penalty from the masses and the stable increment (the "soft constraint" idea: J. O. Hallquist, + *LS-DYNA Theory Manual*, LSTC 2006, the contact chapter's soft-constraint option; the bound is ours: Gershgorin's row sums of the + integrated contact stiffness with the node's directional area, added to the element + frequencies by Weyl's inequality). +- Friction by a regularised (penalty) return map: P. Wriggers, *Computational Contact + Mechanics*, §5.2 and §10.3; J. C. Simo, T. A. Laursen, Comput. Struct. 42(1):97–116, 1992; + rate-dependent coefficients reviewed by J. T. Oden, J. A. C. Martins, CMAME 52:527–634, 1985. + The exponential decay μ = μ_k + (μ_s − μ_k) e^{−d γ̇} is the form the LS-DYNA and Abaqus + manuals document, so the decks carry it unchanged. +- Kinematic rigid walls with Coulomb impulses: T. J. R. Hughes et al., CMAME 8(3):249–276, 1976 + (A5.5), the tangential impulse capped at μ times the normal one. +- Verification: the impact of equal bars, K. F. Graff, *Wave Motion in Elastic Solids*, Dover + 1991, §2.1–2.3, and W. Johnson, *Impact Strength of Materials*, Arnold 1972, §1.6; Hertz's + impact of spheres, H. Hertz, J. reine angew. Math. 92:156–171, 1882, and K. L. Johnson, + *Contact Mechanics*, Cambridge 1985, §11.4; the sliding block, any dynamics text (J. L. Meriam, + L. G. Kraige, *Engineering Mechanics: Dynamics*, §3/5). diff --git a/docs/verification.md b/docs/verification.md index a0f5b0e..ea89b03 100644 --- a/docs/verification.md +++ b/docs/verification.md @@ -59,7 +59,9 @@ cargo xtask oracle --quick # skip cases marked "slow" "known_failure": "optional: why FieldCraft misses the reference today"}]} ``` A measure reads a JSON pointer either from a recipe step's output (`"from": ""`) or from a - command run after the recipe (`"command"`, `"params"`); an optional `"scale"` multiplies it. + command run after the recipe (`"command"`, `"params"`); an optional `"scale"` multiplies it, and + an optional `"over"` (a measure of its own) divides it by another number of the same run, for + ratios such as the stress in one element over the stress in another. A tolerance is `relative` (of the reference), `absolute`, both, or `"exact": true`. A pending case cites its reference values (or `null` until transcribed) and has no recipe yet. - `fixture.json` holds each passing quantity's accepted `value` and its `band` diff --git a/tests/oracle/explicit-bars-collision/README.md b/tests/oracle/explicit-bars-collision/README.md new file mode 100644 index 0000000..2ee19b4 --- /dev/null +++ b/tests/oracle/explicit-bars-collision/README.md @@ -0,0 +1,9 @@ +# Two identical elastic bars colliding (explicit contact) + +Two bars 10 × 0.1 × 0.1 of 200 one-point hex8 elements each (E = 10⁴, ρ = 1, ν = 0, c = 100), +end to end and touching at x = 10 with separate nodes; the first moves at v₀ = 0.1 into the +second, at rest. A frictionless contact pair joins the two end faces (each face's nodes kept out +of the other, segment to segment). One-dimensional theory: the bars stay in contact for 2L/c = +0.2 with a force ρcAv₀/2 across the interface, then the striker stops and the struck bar leaves at +v₀ with all the kinetic energy. Linear bulk viscosity is off (it dissipates at the wave fronts). +Solved to t = 0.3. diff --git a/tests/oracle/explicit-bars-collision/expected.json b/tests/oracle/explicit-bars-collision/expected.json new file mode 100644 index 0000000..b7b2d2f --- /dev/null +++ b/tests/oracle/explicit-bars-collision/expected.json @@ -0,0 +1,107 @@ +{ + "title": "Two identical elastic bars colliding (explicit contact)", + "kind": "theory", + "source": "One-dimensional elastic impact of equal bars: contact for 2L/c, interface force ρcAv₀/2, exchange of velocities, no loss of energy; Graff, Wave Motion in Elastic Solids, Dover 1991, §2.1–2.3; Johnson, Impact Strength of Materials, Arnold 1972, §1.6", + "quantities": [ + { + "id": "duration", + "quantity": "time of the last contact (the bars touch from t = 0) against 2L/c = 0.2", + "measure": { + "from": "solve", + "pointer": "/contacts/0/last_contact" + }, + "value": 0.2, + "tolerance": { + "relative": 0.02 + }, + "reason": "the penalty springs at the interface delay the release by a few element transit times (about 1% of 2L/c with 200 elements per bar)" + }, + { + "id": "force", + "quantity": "the largest force across the interface against ρcAv₀/2 = 0.05", + "measure": { + "from": "solve", + "pointer": "/contacts/0/peak_force" + }, + "value": 0.05, + "tolerance": { + "relative": 0.04 + }, + "reason": "the plateau is ρcAv₀/2 to 0.01%; the peak is the plateau plus the ringing behind the discrete wave front, a few percent" + }, + { + "id": "v_a", + "quantity": "velocity of the struck bar's far end after the impact (0: the striker stops)", + "measure": { + "command": "solve.probe", + "params": { + "step": "Impact", + "at": [ + [ + 0, + 0, + 0 + ] + ] + }, + "pointer": "/nodes/0/velocity/0" + }, + "value": 0.0, + "tolerance": { + "absolute": 0.005 + }, + "reason": "the striker stops on average (its kinetic energy is under 0.01% of the start's, see ke); its far end node keeps the ringing of the discrete release front, a few percent of v₀" + }, + { + "id": "v_b", + "quantity": "velocity of the struck bar's far end after the impact (v₀ = 0.1)", + "measure": { + "command": "solve.probe", + "params": { + "step": "Impact", + "at": [ + [ + 20.0, + 0, + 0 + ] + ] + }, + "pointer": "/nodes/0/velocity/0" + }, + "value": 0.1, + "tolerance": { + "relative": 0.05 + }, + "reason": "the struck bar leaves at v₀ on average; its far end node keeps the ringing of the discrete wave fronts, a few percent" + }, + { + "id": "ke", + "quantity": "kinetic energy after the impact over ½mv₀² (frictionless elastic impact)", + "measure": { + "from": "solve", + "pointer": "/energy/kinetic", + "scale": 2000.0 + }, + "value": 1.0, + "tolerance": { + "relative": 0.01 + }, + "reason": "an elastic impact keeps the kinetic energy; the quadratic bulk viscosity at the compression fronts and the vibration left in the bars take a fraction of a percent" + }, + { + "id": "contact_energy", + "quantity": "energy left in the contact after the separation, over ½mv₀² (its drift)", + "measure": { + "from": "solve", + "pointer": "/energy/contact", + "scale": 2000.0 + }, + "value": 0.0, + "tolerance": { + "absolute": 0.001 + }, + "reason": "the penalty springs give back what they took; the trapezoidal work of a spring that engages between increments drifts by far less than 0.1%" + } + ] +} \ No newline at end of file diff --git a/tests/oracle/explicit-bars-collision/fixture.json b/tests/oracle/explicit-bars-collision/fixture.json new file mode 100644 index 0000000..d32680d --- /dev/null +++ b/tests/oracle/explicit-bars-collision/fixture.json @@ -0,0 +1,29 @@ +{ + "note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.", + "values": { + "contact_energy": { + "band": 0.001, + "value": -0.00007743495029686868 + }, + "duration": { + "band": 2.0242265804764194e-10, + "value": 0.20242265804764192 + }, + "force": { + "band": 5.1469625943477736e-11, + "value": 0.05146962594347773 + }, + "ke": { + "band": 1e-9, + "value": 0.9999532471694228 + }, + "v_a": { + "band": 0.005, + "value": -0.0033987044489733886 + }, + "v_b": { + "band": 1.0339870444897614e-10, + "value": 0.10339870444897613 + } + } +} diff --git a/tests/oracle/explicit-bars-collision/recipe.json b/tests/oracle/explicit-bars-collision/recipe.json new file mode 100644 index 0000000..97fe518 --- /dev/null +++ b/tests/oracle/explicit-bars-collision/recipe.json @@ -0,0 +1,14 @@ +[ + {"command": "mesh.block", "params": {"body": "BarA", "kind": "hex8", "corners": [[0, 0, 0], [10.0, 0, 0], [10.0, 0.1, 0], [0, 0.1, 0], [0, 0, 0.1], [10.0, 0, 0.1], [10.0, 0.1, 0.1], [0, 0.1, 0.1]], "divisions": [200, 1, 1], "faces": {"4": "end"}}}, + {"command": "mesh.block", "params": {"body": "BarB", "kind": "hex8", "corners": [[10.0, 0, 0], [20.0, 0, 0], [20.0, 0.1, 0], [10.0, 0.1, 0], [10.0, 0, 0.1], [20.0, 0, 0.1], [20.0, 0.1, 0.1], [10.0, 0.1, 0.1]], "divisions": [200, 1, 1], "faces": {"6": "start"}, "merge": 0}}, + {"command": "model.set.create", "params": {"name": "A", "geometry": [{"body": "BarA"}]}}, + {"command": "model.set.create", "params": {"name": "B", "geometry": [{"body": "BarB"}]}}, + {"command": "model.set.create", "params": {"name": "a_end", "geometry": [{"body": "BarA", "face": "end"}]}}, + {"command": "model.set.create", "params": {"name": "b_start", "geometry": [{"body": "BarB", "face": "start"}]}}, + {"command": "model.material.create", "params": {"name": "Mat", "elastic": {"type": "isotropic", "young": 10000.0, "poisson": 0.0}, "density": 1.0}}, + {"command": "model.section.create", "params": {"name": "A", "set": "A", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}}, + {"command": "model.section.create", "params": {"name": "B", "set": "B", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}}, + {"command": "interaction.contact.pair", "params": {"name": "impact", "main": "b_start", "secondary": "a_end"}}, + {"command": "model.step.create", "params": {"name": "Impact", "procedure": {"type": "explicit_dynamic", "end_time": 0.3, "initial_velocities": [{"set": "A", "velocity": [0.1, 0, 0]}], "bulk_viscosity": {"linear": 0.0, "quadratic": 1.5}}, "boundary_conditions": [], "loads": []}}, + {"command": "solve.explicit", "params": {"step": "Impact"}, "id": "solve"} +] diff --git a/tests/oracle/explicit-contact-patch/README.md b/tests/oracle/explicit-contact-patch/README.md new file mode 100644 index 0000000..7e1aef4 --- /dev/null +++ b/tests/oracle/explicit-contact-patch/README.md @@ -0,0 +1,9 @@ +# Contact patch test (explicit contact, non-matching meshes) + +Two unit cubes stacked along z, meshed independently: 3 × 3 × 2 hex8 below, 2 × 2 × 2 above, so +the nodes of the interface at z = 1 do not match. A frictionless contact pair joins the lower +block's top and the upper block's bottom. The faces x = 0 and y = 0 of both blocks are symmetry +planes (normal displacement held), the base z = 0 is held in z, and the top z = 2 is pushed down by +0.01 at constant speed over a hundred wave transits (E = 1000, ν = 0.3, ρ = 1). The exact state +is uniaxial stress, uniform in both blocks; the stresses are compared as ratios, since the penalty +takes a small, uniform share of the squeeze. diff --git a/tests/oracle/explicit-contact-patch/expected.json b/tests/oracle/explicit-contact-patch/expected.json new file mode 100644 index 0000000..f302811 --- /dev/null +++ b/tests/oracle/explicit-contact-patch/expected.json @@ -0,0 +1,176 @@ +{ + "title": "Contact patch test: two blocks with non-matching meshes in uniaxial compression (explicit contact)", + "kind": "theory", + "source": "The contact patch test: a contact discretisation must transmit a uniform pressure across a non-matching interface exactly; Taylor & Papadopoulos, On a patch test for contact problems in two dimensions, Nonlinear Computational Mechanics, Springer 1991, 690–702; Crisfield, Re-visiting the contact patch test, IJNME 48:435–449, 2000. Uniaxial stress between frictionless symmetry planes is uniform in both blocks", + "quantities": [ + { + "id": "lower_corner", + "quantity": "σ_zz of the lower block's corner element under the interface (element 10) over that of the upper block's corner element over it (element 19)", + "measure": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 10 + ] + }, + "pointer": "/elements/0/stress/2", + "over": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 19 + ] + }, + "pointer": "/elements/0/stress/2" + } + }, + "value": 1.0, + "tolerance": { + "relative": 0.005 + }, + "reason": "a uniform stress is an equilibrium of the discrete system (the integrated contact forces are the consistent nodal forces of a uniform pressure on both meshes); what is left is the ringing of the loading, a hundred wave transits slow" + }, + { + "id": "lower_centre", + "quantity": "σ_zz of the lower block's centre element under the interface (element 14) over that of element 19", + "measure": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 14 + ] + }, + "pointer": "/elements/0/stress/2", + "over": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 19 + ] + }, + "pointer": "/elements/0/stress/2" + } + }, + "value": 1.0, + "tolerance": { + "relative": 0.005 + }, + "reason": "as for the corner element" + }, + { + "id": "lower_base", + "quantity": "σ_zz of the lower block's bottom corner element (element 1) over that of element 19", + "measure": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 1 + ] + }, + "pointer": "/elements/0/stress/2", + "over": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 19 + ] + }, + "pointer": "/elements/0/stress/2" + } + }, + "value": 1.0, + "tolerance": { + "relative": 0.005 + }, + "reason": "as for the corner element" + }, + { + "id": "upper_top", + "quantity": "σ_zz of the upper block's top far corner element (element 26) over that of element 19", + "measure": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 26 + ] + }, + "pointer": "/elements/0/stress/2", + "over": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 19 + ] + }, + "pointer": "/elements/0/stress/2" + } + }, + "value": 1.0, + "tolerance": { + "relative": 0.005 + }, + "reason": "as for the corner element" + }, + { + "id": "force", + "quantity": "the force across the interface over −σ_zz A (A = 1)", + "measure": { + "from": "solve", + "pointer": "/contacts/0/peak_force", + "over": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 19 + ] + }, + "pointer": "/elements/0/stress/2", + "scale": -1.0 + } + }, + "value": 1.0, + "tolerance": { + "relative": 0.01 + }, + "reason": "the integrated contact pressure carries the whole load; the largest force of the run is the final one in a slow squeeze, up to the ringing" + }, + { + "id": "lateral", + "quantity": "σ_xx of element 14 over its σ_zz (uniaxial stress: 0)", + "measure": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 14 + ] + }, + "pointer": "/elements/0/stress/0", + "over": { + "command": "solve.probe", + "params": { + "step": "Squeeze", + "elements": [ + 14 + ] + }, + "pointer": "/elements/0/stress/2" + } + }, + "value": 0.0, + "tolerance": { + "absolute": 0.005 + }, + "reason": "frictionless symmetry planes and a frictionless interface leave the blocks free to expand sideways; the non-matching interface must not hold them" + } + ] +} \ No newline at end of file diff --git a/tests/oracle/explicit-contact-patch/fixture.json b/tests/oracle/explicit-contact-patch/fixture.json new file mode 100644 index 0000000..41fd27b --- /dev/null +++ b/tests/oracle/explicit-contact-patch/fixture.json @@ -0,0 +1,29 @@ +{ + "note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.", + "values": { + "force": { + "band": 1.0025883998151134e-9, + "value": 1.0025883998151133 + }, + "lateral": { + "band": 0.005, + "value": -0.00036052361311792505 + }, + "lower_base": { + "band": 1.00013429984765e-9, + "value": 1.00013429984765 + }, + "lower_centre": { + "band": 1.0002867148178652e-9, + "value": 1.0002867148178651 + }, + "lower_corner": { + "band": 1.0001179568787019e-9, + "value": 1.0001179568787018 + }, + "upper_top": { + "band": 1e-9, + "value": 0.9994358602920412 + } + } +} diff --git a/tests/oracle/explicit-contact-patch/recipe.json b/tests/oracle/explicit-contact-patch/recipe.json new file mode 100644 index 0000000..b4f3ff7 --- /dev/null +++ b/tests/oracle/explicit-contact-patch/recipe.json @@ -0,0 +1,22 @@ +[ + {"command": "mesh.block", "params": {"body": "Lower", "kind": "hex8", "corners": [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]], "divisions": [3, 3, 2], "faces": {"1": "bottom", "2": "top", "3": "y0", "6": "x0"}}}, + {"command": "mesh.block", "params": {"body": "Upper", "kind": "hex8", "corners": [[0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1], [0, 0, 2], [1, 0, 2], [1, 1, 2], [0, 1, 2]], "divisions": [2, 2, 2], "faces": {"1": "bottom", "2": "top", "3": "y0", "6": "x0"}, "merge": 0}}, + {"command": "model.set.create", "params": {"name": "lower", "geometry": [{"body": "Lower"}]}}, + {"command": "model.set.create", "params": {"name": "upper", "geometry": [{"body": "Upper"}]}}, + {"command": "model.set.create", "params": {"name": "lower_top", "geometry": [{"body": "Lower", "face": "top"}]}}, + {"command": "model.set.create", "params": {"name": "upper_bottom", "geometry": [{"body": "Upper", "face": "bottom"}]}}, + {"command": "model.set.create", "params": {"name": "x0", "geometry": [{"body": "Lower", "face": "x0"}, {"body": "Upper", "face": "x0"}]}}, + {"command": "model.set.create", "params": {"name": "y0", "geometry": [{"body": "Lower", "face": "y0"}, {"body": "Upper", "face": "y0"}]}}, + {"command": "model.set.create", "params": {"name": "z0", "geometry": [{"body": "Lower", "face": "bottom"}]}}, + {"command": "model.set.create", "params": {"name": "top", "geometry": [{"body": "Upper", "face": "top"}]}}, + {"command": "model.material.create", "params": {"name": "Mat", "elastic": {"type": "isotropic", "young": 1000.0, "poisson": 0.3}, "density": 1.0}}, + {"command": "model.section.create", "params": {"name": "lower", "set": "lower", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}}, + {"command": "model.section.create", "params": {"name": "upper", "set": "upper", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}}, + {"command": "interaction.contact.pair", "params": {"name": "interface", "main": "lower_top", "secondary": "upper_bottom"}}, + {"command": "model.boundary_condition.create", "params": {"name": "x0", "type": "displacement", "target": {"set": "x0"}, "values": [0, null, null, null, null, null]}}, + {"command": "model.boundary_condition.create", "params": {"name": "y0", "type": "displacement", "target": {"set": "y0"}, "values": [null, 0, null, null, null, null]}}, + {"command": "model.boundary_condition.create", "params": {"name": "z0", "type": "displacement", "target": {"set": "z0"}, "values": [null, null, 0, null, null, null]}}, + {"command": "model.boundary_condition.create", "params": {"name": "push", "type": "displacement", "target": {"set": "top"}, "values": [null, null, -0.01, null, null, null]}}, + {"command": "model.step.create", "params": {"name": "Squeeze", "procedure": {"type": "explicit_dynamic", "end_time": 5.451}, "boundary_conditions": ["x0", "y0", "z0", "push"], "loads": []}}, + {"command": "solve.explicit", "params": {"step": "Squeeze"}, "id": "solve"} +] diff --git a/tests/oracle/explicit-hertz-spheres/README.md b/tests/oracle/explicit-hertz-spheres/README.md new file mode 100644 index 0000000..6071145 --- /dev/null +++ b/tests/oracle/explicit-hertz-spheres/README.md @@ -0,0 +1,8 @@ +# Two elastic spheres colliding head-on (explicit contact, Hertz) + +Two identical steel spheres, R = 10 mm, E = 210 GPa, ν = 0.3, ρ = 7.85 t/m³ (mm, t, s, N, MPa), +meet head-on at a relative speed of 1 m/s. Hertz's quasi-static theory gives the approach +δ = 0.020441 mm, the contact radius a = √(R* δ) = 0.320 mm, the peak force 1005.4 N and the +contact duration 60.16 μs. The case waits for a graded sphere mesh built from commands: the contact +radius is a thirtieth of the sphere's, so the mesh must be fine at the poles and coarse elsewhere. +The explicit crate's tests check the same theory on a cubed-sphere mesh generated in code. diff --git a/tests/oracle/explicit-hertz-spheres/expected.json b/tests/oracle/explicit-hertz-spheres/expected.json new file mode 100644 index 0000000..2eecb01 --- /dev/null +++ b/tests/oracle/explicit-hertz-spheres/expected.json @@ -0,0 +1,11 @@ +{ + "title": "Two elastic spheres in a head-on collision: Hertz's force and contact duration (explicit contact)", + "kind": "theory", + "source": "Hertz's quasi-static impact: with R* = R/2, m* = m/2, E* = E/(2(1 − ν²)), the approach reaches δ = (15 m* v²/(16 E* √R*))^(2/5), the force (4/3) E* √R* δ^(3/2), and the contact lasts 2.9432 δ/v; Hertz, J. reine angew. Math. 92:156–171, 1882; K. L. Johnson, Contact Mechanics, Cambridge 1985, §11.4", + "status": "pending", + "requires": "a graded sphere mesh from commands (a cubed-sphere hex mesher, or tet refinement towards a point); until then the explicit crate's tests `spheres_collide_as_hertz_says` (a graded cubed-sphere quarter of each sphere, 27 648 hex8: peak force +2.3%, duration −0.4%) and its coarse twin in the main gate check this case on a mesh generated in code", + "quantities": [ + {"id": "force", "quantity": "the largest contact force, two steel spheres R = 10 mm (E = 210 GPa, ν = 0.3, ρ = 7.85 t/m³) meeting at 1 m/s (N)", "value": 1005.38}, + {"id": "duration", "quantity": "the contact duration (s)", "value": 6.0163e-5} + ] +} diff --git a/tests/oracle/explicit-sliding-block-contact/README.md b/tests/oracle/explicit-sliding-block-contact/README.md new file mode 100644 index 0000000..41449aa --- /dev/null +++ b/tests/oracle/explicit-sliding-block-contact/README.md @@ -0,0 +1,8 @@ +# Block sliding to a stop on a deformable plate (explicit contact, decaying friction) + +The block of `explicit-sliding-block-wall` (1 × 1 × 0.5, 4 × 4 × 2 hex8, E = 10⁵, ν = 0.3, +ρ = 1) rests on a plate 4.5 × 2 × 0.5 of 9 × 4 × 1 hex8 held at its base, the two meshes not +matching; a contact pair joins the plate's top and the block's bottom with Coulomb friction whose +coefficient falls with the slip rate, μ = 0.2 + 0.3 e^(−1.5 v). Under gravity g = 10 the block +is launched along x at v₀ = 3. A rigid block decelerates at g μ(v), stops after 1.2001 and travels +∫₀^v₀ v dv / (g μ(v)) = 2.01524. Solved to t = 1.56. diff --git a/tests/oracle/explicit-sliding-block-contact/expected.json b/tests/oracle/explicit-sliding-block-contact/expected.json new file mode 100644 index 0000000..8e60ccd --- /dev/null +++ b/tests/oracle/explicit-sliding-block-contact/expected.json @@ -0,0 +1,15 @@ +{ + "title": "Block sliding to a stop on a deformable plate, friction decaying with the slip rate (explicit contact)", + "kind": "theory", + "source": "A rigid block launched at v₀ with friction μ(v) = μk + (μs − μk) e^(−d v) decelerates at g μ(v); it travels ∫₀^v₀ v dv / (g μ(v)) (Simpson's rule, 200 000 intervals) and stops after [v₀ + ln((μk + (μs − μk) e^(−d v₀))/μs)/d]/(g μk) = 1.2001; Meriam & Kraige, Engineering Mechanics: Dynamics, §3/5; the exponential decay law as in Oden & Martins, CMAME 52:527–634, 1985", + "quantities": [ + {"id": "distance", "quantity": "distance travelled (top centre node) against ∫ v dv / (g μ(v)) = 2.01524", + "measure": {"command": "solve.probe", "params": {"step": "Slide", "at": [[0.5, 0.5, 0.5]]}, "pointer": "/nodes/0/displacement/0"}, + "value": 2.015240387808855, "tolerance": {"relative": 0.01}, + "reason": "the block bounces on the penalty springs and the plate under it shears elastically, so the normal force and the slip rate it feels vary by a little around their rigid-body values"}, + {"id": "velocity", "quantity": "velocity of the top centre node after the stop", + "measure": {"command": "solve.probe", "params": {"step": "Slide", "at": [[0.5, 0.5, 0.5]]}, "pointer": "/nodes/0/velocity/0"}, + "value": 0.0, "tolerance": {"absolute": 0.03}, + "reason": "the block sticks; what is left is its elastic vibration, under 1% of v₀"} + ] +} diff --git a/tests/oracle/explicit-sliding-block-contact/fixture.json b/tests/oracle/explicit-sliding-block-contact/fixture.json new file mode 100644 index 0000000..3935c3a --- /dev/null +++ b/tests/oracle/explicit-sliding-block-contact/fixture.json @@ -0,0 +1,13 @@ +{ + "note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.", + "values": { + "distance": { + "band": 2.0152403878088553e-9, + "value": 2.011776866282693 + }, + "velocity": { + "band": 0.03, + "value": -0.00003285207146964345 + } + } +} diff --git a/tests/oracle/explicit-sliding-block-contact/recipe.json b/tests/oracle/explicit-sliding-block-contact/recipe.json new file mode 100644 index 0000000..d21aafc --- /dev/null +++ b/tests/oracle/explicit-sliding-block-contact/recipe.json @@ -0,0 +1,17 @@ +[ + {"command": "mesh.block", "params": {"body": "Block", "kind": "hex8", "corners": [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 0.5], [1, 0, 0.5], [1, 1, 0.5], [0, 1, 0.5]], "divisions": [4, 4, 2], "faces": {"1": "bottom"}}}, + {"command": "mesh.block", "params": {"body": "Plate", "kind": "hex8", "corners": [[-0.5, -0.5, -0.5], [4.0, -0.5, -0.5], [4.0, 1.5, -0.5], [-0.5, 1.5, -0.5], [-0.5, -0.5, 0], [4.0, -0.5, 0], [4.0, 1.5, 0], [-0.5, 1.5, 0]], "divisions": [9, 4, 1], "faces": {"1": "base", "2": "top"}, "merge": 0}}, + {"command": "model.set.create", "params": {"name": "block", "geometry": [{"body": "Block"}]}}, + {"command": "model.set.create", "params": {"name": "plate", "geometry": [{"body": "Plate"}]}}, + {"command": "model.set.create", "params": {"name": "bottom", "geometry": [{"body": "Block", "face": "bottom"}]}}, + {"command": "model.set.create", "params": {"name": "top", "geometry": [{"body": "Plate", "face": "top"}]}}, + {"command": "model.set.create", "params": {"name": "base", "geometry": [{"body": "Plate", "face": "base"}]}}, + {"command": "model.material.create", "params": {"name": "Mat", "elastic": {"type": "isotropic", "young": 100000.0, "poisson": 0.3}, "density": 1.0}}, + {"command": "model.section.create", "params": {"name": "block", "set": "block", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}}, + {"command": "model.section.create", "params": {"name": "plate", "set": "plate", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}}, + {"command": "interaction.contact.pair", "params": {"name": "slide", "main": "top", "secondary": "bottom", "friction": {"static": 0.5, "kinetic": 0.2, "decay": 1.5}}}, + {"command": "model.boundary_condition.create", "params": {"name": "hold", "type": "fixed", "target": {"set": "base"}}}, + {"command": "model.load.create", "params": {"name": "weight", "type": "gravity", "acceleration": [0, 0, -10.0], "set": "block"}}, + {"command": "model.step.create", "params": {"name": "Slide", "procedure": {"type": "explicit_dynamic", "end_time": 1.56, "initial_velocities": [{"set": "block", "velocity": [3.0, 0, 0]}]}, "boundary_conditions": ["hold"], "loads": ["weight"]}}, + {"command": "solve.explicit", "params": {"step": "Slide"}, "id": "solve"} +] diff --git a/tests/oracle/explicit-sliding-block-wall/README.md b/tests/oracle/explicit-sliding-block-wall/README.md new file mode 100644 index 0000000..5d7de46 --- /dev/null +++ b/tests/oracle/explicit-sliding-block-wall/README.md @@ -0,0 +1,7 @@ +# Block sliding to a stop on a rigid wall (explicit, Coulomb friction) + +A block 1 × 1 × 0.5 of 4 × 4 × 2 one-point hex8 elements (E = 10⁵, ν = 0.3, ρ = 1, so m = 0.5) +rests on a rigid wall at z = 0 (its bottom nodes) with Coulomb friction μ = 0.3, under gravity +g = 10, and is launched along x at v₀ = 3. Theory for a rigid block: deceleration μg, stop after +v₀/(μg) = 1 at v₀²/(2μg) = 1.5, the kinetic energy ½mv₀² = 2.25 all dissipated by friction. +Solved to t = 1.3. diff --git a/tests/oracle/explicit-sliding-block-wall/expected.json b/tests/oracle/explicit-sliding-block-wall/expected.json new file mode 100644 index 0000000..ad64576 --- /dev/null +++ b/tests/oracle/explicit-sliding-block-wall/expected.json @@ -0,0 +1,19 @@ +{ + "title": "Block sliding to a stop on a rigid wall with Coulomb friction (explicit)", + "kind": "theory", + "source": "A rigid block launched at v₀ on a horizontal plane with Coulomb friction μ decelerates at μg and stops after v₀/(μg) at v₀²/(2μg), its kinetic energy all dissipated by friction; Meriam & Kraige, Engineering Mechanics: Dynamics, §3/5 (work and energy)", + "quantities": [ + {"id": "distance", "quantity": "distance travelled (top centre node) against v₀²/(2μg) = 1.5", + "measure": {"command": "solve.probe", "params": {"step": "Slide", "at": [[0.5, 0.5, 0.5]]}, "pointer": "/nodes/0/displacement/0"}, + "value": 1.5, "tolerance": {"relative": 0.01}, + "reason": "the wall's friction impulses are exact Coulomb impulses; the block's elastic vibration from the sudden gravity and the friction moment move the node by a fraction of a percent"}, + {"id": "velocity", "quantity": "velocity of the top centre node after the stop", + "measure": {"command": "solve.probe", "params": {"step": "Slide", "at": [[0.5, 0.5, 0.5]]}, "pointer": "/nodes/0/velocity/0"}, + "value": 0.0, "tolerance": {"absolute": 0.03}, + "reason": "the block sticks; what is left is its elastic vibration, under 1% of v₀"}, + {"id": "dissipated", "quantity": "friction dissipation over ½mv₀² = 2.25", + "measure": {"from": "solve", "pointer": "/energy/friction", "scale": 0.4444444444444444}, + "value": 1.0, "tolerance": {"relative": 0.01}, + "reason": "the kinetic energy goes into friction; the vibration and the work of gravity as the block settles add or take a fraction of a percent"} + ] +} diff --git a/tests/oracle/explicit-sliding-block-wall/fixture.json b/tests/oracle/explicit-sliding-block-wall/fixture.json new file mode 100644 index 0000000..4e321ea --- /dev/null +++ b/tests/oracle/explicit-sliding-block-wall/fixture.json @@ -0,0 +1,17 @@ +{ + "note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.", + "values": { + "dissipated": { + "band": 1.0012078586619573e-9, + "value": 1.0012078586619573 + }, + "distance": { + "band": 1.5000000000000002e-9, + "value": 1.4999205464422256 + }, + "velocity": { + "band": 0.03, + "value": 0.00027058229054461555 + } + } +} diff --git a/tests/oracle/explicit-sliding-block-wall/recipe.json b/tests/oracle/explicit-sliding-block-wall/recipe.json new file mode 100644 index 0000000..51fc7a1 --- /dev/null +++ b/tests/oracle/explicit-sliding-block-wall/recipe.json @@ -0,0 +1,10 @@ +[ + {"command": "mesh.block", "params": {"body": "Block", "kind": "hex8", "corners": [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 0.5], [1, 0, 0.5], [1, 1, 0.5], [0, 1, 0.5]], "divisions": [4, 4, 2], "faces": {"1": "bottom"}}}, + {"command": "model.set.create", "params": {"name": "block", "geometry": [{"body": "Block"}]}}, + {"command": "model.set.create", "params": {"name": "bottom", "geometry": [{"body": "Block", "face": "bottom"}]}}, + {"command": "model.material.create", "params": {"name": "Mat", "elastic": {"type": "isotropic", "young": 100000.0, "poisson": 0.3}, "density": 1.0}}, + {"command": "model.section.create", "params": {"name": "block", "set": "block", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}}, + {"command": "model.load.create", "params": {"name": "weight", "type": "gravity", "acceleration": [0, 0, -10.0], "set": "block"}}, + {"command": "model.step.create", "params": {"name": "Slide", "procedure": {"type": "explicit_dynamic", "end_time": 1.3, "initial_velocities": [{"set": "block", "velocity": [3.0, 0, 0]}], "rigid_walls": [{"set": "bottom", "point": [0, 0, 0], "normal": [0, 0, 1], "friction": {"static": 0.3}}]}, "boundary_conditions": [], "loads": ["weight"]}}, + {"command": "solve.explicit", "params": {"step": "Slide"}, "id": "solve"} +]