From 15a9ac9157358bab3812c8c99b63285dfdb958c4 Mon Sep 17 00:00:00 2001 From: Irakli Kourtis Date: Tue, 13 Oct 2026 21:36:46 -0400 Subject: [PATCH] 1D elements I: trusses (rod2, rod3), point masses, rotary inertia, springs and dashpots (#38) * Model, fem and decks: point masses, rotary inertia, springs, dashpots and rod3 trusses Model: PointMass and Spring definitions (spring | dashpot; to ground or two-point; along a degree of freedom, a vector or the line of the ends) on a set or reference point, resolved per node; ElementKind::Rod3. fem: rod2/rod3 truss kernel on Basis::Line2/Line3 declaring ElementDofs (translations). Decks: Abaqus T3D3, MASS, ROTARYI, SPRING1/2/A, DASHPOT1/2/A; Nastran CONM1, CELAS2, CDAMP2; LS-DYNA ELEMENT_MASS/INERTIA and ELEMENT_DISCRETE; what a format cannot take is reported. Commands model.mass.* and model.spring.*. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * Implicit solver: trusses, point masses, rotary inertia and springs on the DofLayout Trusses declare their translations like solids; point masses, inertia and springs ask the layout for the kinds they need at their nodes (translations, rotations), and reference points they use get six degrees of freedom. The matrix pattern takes rows for those points; dashpots are ignored with a warning. Existing results are bit-identical (hash-checked). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * Explicit solver: trusses, point masses, springs and dashpots on the node-major layout Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * docs: fem.md on trusses and discrete elements in the layout; regenerate features.md Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * docs: trusses in plane models and corotational explicit trusses Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm --------- Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com> Co-authored-by: Claude Opus 5.5 --- ATTRIBUTION.md | 4 +- ROADMAP.md | 4 +- assets/icons/interaction.mass.vectorcraft | Bin 1886 -> 0 bytes assets/icons/interaction.spring.vectorcraft | Bin 1897 -> 0 bytes ...raction.mass.svg => model.mass.create.svg} | 2 +- assets/icons/model.mass.create.vectorcraft | Bin 0 -> 1886 bytes ...ion.spring.svg => model.spring.create.svg} | 2 +- assets/icons/model.spring.create.vectorcraft | Bin 0 -> 1896 bytes ...ction.mass.json => model.mass.create.json} | 0 ...n.spring.json => model.spring.create.json} | 0 crates/decks/src/abaqus.rs | 149 +++- crates/decks/src/lib.rs | 27 +- crates/decks/src/lsdyna/elements.rs | 4 +- crates/decks/src/lsdyna/mod.rs | 138 +++- crates/decks/src/nastran/mod.rs | 62 +- crates/decks/tests/discrete.rs | 166 +++++ crates/engine/src/catalogue.rs | 18 +- crates/engine/src/cmd/define.rs | 35 +- crates/engine/src/cmd/mod.rs | 6 + crates/engine/src/cmd/solve.rs | 28 + crates/engine/src/lib.rs | 7 + crates/explicit/src/discrete.rs | 316 +++++++++ crates/explicit/src/elements.rs | 16 + crates/explicit/src/integrate.rs | 40 +- crates/explicit/src/lib.rs | 1 + crates/explicit/src/system.rs | 90 ++- crates/explicit/tests/one_d.rs | 307 +++++++++ crates/fem/src/basis.rs | 1 + crates/fem/src/lib.rs | 3 + crates/fem/src/truss.rs | 404 +++++++++++ crates/io/src/mesh.rs | 3 +- crates/mesh/src/quality/mod.rs | 4 +- crates/model/src/check.rs | 151 +++- crates/model/src/constraints.rs | 1 + crates/model/src/discrete.rs | 501 ++++++++++++++ crates/model/src/lib.rs | 25 +- crates/model/src/mesh.rs | 20 +- crates/model/src/sets.rs | 1 + crates/render/src/display.rs | 2 +- crates/solver/src/analysis.rs | 85 ++- crates/solver/src/constraints.rs | 7 +- crates/solver/src/discrete.rs | 198 ++++++ crates/solver/src/lib.rs | 1 + crates/solver/src/sparse.rs | 14 +- crates/solver/src/system.rs | 176 ++++- crates/solver/tests/one_d.rs | 652 ++++++++++++++++++ crates/ui-egui/src/theme.rs | 2 +- docs/architecture.md | 26 +- docs/elements.md | 12 +- docs/features.md | 6 +- docs/fem.md | 6 + docs/parity.md | 3 +- docs/solvers/literature.md | 10 + 53 files changed, 3626 insertions(+), 110 deletions(-) delete mode 100644 assets/icons/interaction.mass.vectorcraft delete mode 100644 assets/icons/interaction.spring.vectorcraft rename assets/icons/{interaction.mass.svg => model.mass.create.svg} (95%) create mode 100644 assets/icons/model.mass.create.vectorcraft rename assets/icons/{interaction.spring.svg => model.spring.create.svg} (95%) create mode 100644 assets/icons/model.spring.create.vectorcraft rename assets/icons/src/{interaction.mass.json => model.mass.create.json} (100%) rename assets/icons/src/{interaction.spring.json => model.spring.create.json} (100%) create mode 100644 crates/decks/tests/discrete.rs create mode 100644 crates/explicit/src/discrete.rs create mode 100644 crates/explicit/tests/one_d.rs create mode 100644 crates/fem/src/truss.rs create mode 100644 crates/model/src/discrete.rs create mode 100644 crates/solver/src/discrete.rs create mode 100644 crates/solver/tests/one_d.rs diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index 6912350..de90298 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -39,9 +39,9 @@ is committed and listed here with its exports. | assets/icons/interaction.contact.general.svg, assets/icons/interaction.contact.general.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.contact.general.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/interaction.contact.pair.svg, assets/icons/interaction.contact.pair.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.contact.pair.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/interaction.fsi.svg, assets/icons/interaction.fsi.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.fsi.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | -| assets/icons/interaction.mass.svg, assets/icons/interaction.mass.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.mass.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | +| assets/icons/model.mass.create.svg, assets/icons/model.mass.create.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.mass.create.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/interaction.rigid_wall.svg, assets/icons/interaction.rigid_wall.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.rigid_wall.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | -| assets/icons/interaction.spring.svg, assets/icons/interaction.spring.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.spring.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | +| assets/icons/model.spring.create.svg, assets/icons/model.spring.create.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.spring.create.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/interaction.tie.svg, assets/icons/interaction.tie.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/interaction.tie.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/mesh.delaunay.points.svg, assets/icons/mesh.delaunay.points.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.delaunay.points.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/mesh.element.add.svg, assets/icons/mesh.element.add.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/mesh.element.add.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | diff --git a/ROADMAP.md b/ROADMAP.md index 4cb66c1..da4cbb7 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -22,10 +22,10 @@ hex meshes today and it gives the best elements where stresses are read. | M1 Geometry in and surface preparation | Open `.solidcraft` designs with persistent face names, assemblies placed (done); face-tagged 3MF (done); test corpus, `cargo xtask corpus` (done); exact predicates (done); quality metrics and certified validity gate (done); weld, closed/manifold/orientation and self-intersection checks, feature edges (`modelcraft-mesh::surface`, `mesh.surface.check`, done) | in progress | | M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity; the `modelcraft_mesh::size::SizeField` trait and edge-division helpers are in); Bowyer–Watson Delaunay (exact, symbolic perturbation, verified; `modelcraft-mesh::delaunay`, done); robust CDT boundary recovery (Diazzi et al. 2023: implicit LNC Steiner points on segments, gift-wrapping facet recovery, inside/outside labelling, B-rep face tags on boundary tet faces, few-ulp rounding repair checked by the validity gate; `modelcraft-mesh::cdt` and `mesh.tet.generate`, done: all 20 valid corpus parts); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress | | M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input (done: `deck.write`), Nastran SOL 101/103 (done: `deck.write`) and LS-DYNA implicit (done); the `.modelcraft` file (text definition, binary mesh, geometry references reopened on load, journal; atomic saves, autosave and recovery; undo/redo on the journal: done) | in progress | -| M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion; 2D elements: plane stress, plane strain and axisymmetric tri3, tri6, quad4 (incompatible modes QM6), quad8 with two degrees of freedom per node, centrifugal loads, `mesh.planar` 2D meshes of polygons and planar faces, verified by patch tests, Kirsch, Lamé (plane strain and axisymmetric), rotating discs and cylinders, and a plane stress cantilever; deck mapping to Abaqus CPS/CPE/CAX, Nastran membranes, LS-DYNA ELFORM 12/13); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started | +| M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion; 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); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started | | M5 Structured hex | 2D quad meshing of faces (`modelcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | | M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned | -| M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, springs, point masses; amplitudes, multiple steps; nonlinear static steps in decks | planned | +| M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, connectors; amplitudes, multiple steps; nonlinear static steps in decks | 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; heat transfer | planned | | 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 — `modelcraft-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 | diff --git a/assets/icons/interaction.mass.vectorcraft b/assets/icons/interaction.mass.vectorcraft deleted file mode 100644 index 47cfda1e0499bd2de31fd9dff7dbc53ceb5e4921..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1886 zcmV-k2ch^MiwFP!00000|K(ciZre5#eicK%4MUWuTY&!Px}q!Etwl1l7=pl{ZN}yz zi<&M;yuk0i!<%e5PDatmHedk)OC-hf={px59_qon=GB@s-kEnx=MAss6}fJ_G2B%( zr)Gzj5#CT?+<1t?0k+wHiMRW}vCp^nnhSo|egTAIc|f+H;~(@;Q| zcgGgZ3do*7=bV(nFnC*t*?;G(=x(=H+ehq|1fCqLA84)HnYdJR0f8&7l?jtxMqnQS zn)p}Z{X|0o4cXp3deo40kAQ8=b6(ShQXF8BGzpD`t{N+NLA`&E?c{ao&NP{0kGGA! zPbMJ_5ss3`kNhMKc1#D*8_UrG!X@(iBT=S%)m$i?#8I3^IKUasf-u#+n#YSZ=YX)A zJnEpR$)m?mP#X1=tRuos%lI|@GW4_553^7!FGf;dMtQ2!&M%cs^|XAdB9REOg0ns4 zKtnu1K8|AGKZ|1&rSVQXvOk=pS%e|R_n4R>ItVijEJ z-;zQ9rBFv0`*QmGYatyPeEl!3jlz6k@9+=f>O0cRS5#SnUIQN#CF%FKccLliLK z9H|%(Xaoq$2$5O?`2t~UNk=E(!I6M~P3jN{d?&(jg%X>!O+txGPNSggt#npU-mYgd ze{({T!V(f19mScCvKVO~f#?aM$QsCR;KuyU83`nEky$^oQ2t@ViMFMIP(~jr3W=#> zHqtYfG!aZW*B}C0Qz0-p$Pq(66L97=tqdSWR6fR|r%DmaiY%BoUK#H$&$^1=XI4P* z5&G(em#QoSp+TKAh*f7;*#zoPC7H8kqnp=kQE~|f%t1glbED_yK2{uger{2~g3Z4; z0a1hPtQ+_s`RF7D)g1p<1~1y8KAQ8rmf2aBM_!Ic%RJ~c$o(>>M-zqD(U?dI5W)jZ|3s@M5^ z0(FZhLaIj8C5$r1js|QxC6E(}@jEPX!SmC>r$!pXOVkqk5U{E&4h3>b8Lhue=@P*3 z(>n#^q`FnwE5Dm$E3-N-lAi{E<(DfqpVrML?_2MW^qRD}(3l)Cx)3A0*ZW}L=t1a9 zeMfTiw6pu1z?X}-{x(W{-8G%@dwn*2aN0Wh(l$AyINpEb`xx8z3Z8kdpD4NAj6TW& z*(-jfWnJ(6C%D0T+=S0?BQ_YTQ`G%(C$29=i}j(mrYhK@3O}MsWVD<0jOT(35l>pS zpys+_^7EWpWX}aiOT-r?;=7y$yJg}_OBk5K)D@}4_ly*(g7c=HmL)S|*&6EDFzNwD Y74*KDiWx;lzU$cj1-hG4PX{Uh01s}Kga7~l diff --git a/assets/icons/interaction.spring.vectorcraft b/assets/icons/interaction.spring.vectorcraft deleted file mode 100644 index caf2653700d3422f798bf1cd7f63d6fae5caa7ce..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1897 zcmV-v2bTCBiwFP!00000|K(b3Z`(Ey{wsz)3x-6!2iS)$E4rfHS|me@p$H6GqHQj+ zsOcq%7x=&LNQzFRC@w`O!+;uaVv!W@(cN>$<5Bn4id9?OSZCG^NgGzB6<#%#4Odl7 zSz(>|wzVb&sW5cHs4rsFIUhB1=E>ZFu^CI-Eh%7<`-l5n+iGZ&6S#*K4cwnL&_1ur ziWY0QsV!(zL+59DNAemPQkJu7v1J**}4CT1=0mXCPT_03!u zSOF?7{m!w@oVjCLJN4tvwa$<)nz~-2Vs&S8@&=mC$m< zRTYQ)6M?bCRUD$o4{U46uM#&3{m2RT`@MaTRW}v8CL_`cRuBuG@RnxH1{4JvzuG@h zPKlw@mj9r&ZkPP0d>VMZ4eNAW!e4ieY!}TPumUX}u9L)fT@?69m;^owp_jFd>ISh+ z5XbNkj_3G6bfpkjms`d_D1Q__xd+^EojFpIqrr`?ue#;z-au}wy2q#Y=W!3axE=u~ zC+|i-eX=9iy?K=MK1$RG(y_5p^4MdKjUE7jmNaKI$s}uqrE_{C=vi;DFg3|@TGpH@ zx7L|A5BhV4wzU~=E-BpojeD^dIzoeYD`8&8cDk=CzDSxt(vy3N_rrQhGIW*XJHC=k zau(zcy~P>a+IJQh-;?YecNT;O-Xa^Gbr&^W*Iy+3XB}n)?GcB0`sA$sNZB(dWqUguk(E8I!>mh9X=b{O!x6qh_-U~R!?|BG=;c9*Am zlL)Vzjfy|vupt5j(W7fnYhb&*7H>Etd@>{+PAhpTS{@OmOGUtH@P+pJ7B}gJNW5J# z^6hsp$kmZpv(DKP=MM=_$r?4YeVFCz2JCNSLk^r9PyAoXn5HZB=NzZkhRks~$~YZa z^rc>q%ur^6c(91nO5jR!QIHpiSwJZuK zfm{_8l*1ya3FXapDH3TTB(d+gp0ZJxIBpWU9T3EOP!MPX`4d7jyER4~E%C(1 zS32pqG!e?|T}5DrV3^4C=n_Xbi+Ca{CxH!pWltw}k-BcV2g+jbZFiiCk3R@tBqhmt{q}z<>8|I(3BBzo= z#8;&JHSz>sCoKA{DFRn4w{ED`-$W2orHjgMw2{0oGP-7_ScDzCKoM1Xj29`(cpW%x zG+!voNn?mcp>TVlRCsNquaM7eRp!P)U5S=JL|Ap>0}`PheaHO`Lk+Sq4I&|gZ4ej@ zazu|K4D_jb!D>@^(epJRti)$}bx3)4h-n6&mO;zuH+>U(d+ZcRg3B;EHRH51pK zimCbwoN+K!<_JWjK!-&|CTjZOoUK;Glo^AZUFC?97E^Qf^9=)8xG~;;D@Ckjs?l5Z zWtyfB<`NgCxlhy!9{nXpmUvXZv#;Z%v+sUB5<4)dazJ((pQrD&{uh-)WlA&Ia){$J9pcr2HuxHPeBGwho zC_m}4AAX)4YF{5|0VOo(s~c9zvJ8XqIFy(j!J>~<%Xt - interaction.mass + model.mass.create diff --git a/assets/icons/model.mass.create.vectorcraft b/assets/icons/model.mass.create.vectorcraft new file mode 100644 index 0000000000000000000000000000000000000000..c92d8163a6f6c71f9aa94a1576bb1e42715b16d4 GIT binary patch literal 1886 zcmV-k2ch^MiwFP!00000|K(b3Z`(Ey{wsz)8-^%RZvpzD>x!;uw-(9JVhDnOmT8-d zENXg5;syTiJHFYL<75<_Yy%b`utZY4htD04$0HxSIj>iw_0GIonzy{pYcg-WF}&3c z6% z*j$ahf~;wM49I6uAyM40awbr;>q=C}FU-J1`PZ5@6mSX%EnQWGV8V#5YF4&nD)?&a z4hr%mD`qtT{Twnl6wCh+sW(FooO=19`9QH zo=id@J z^Joi2bsjy2g3_p`WbF}lTE?$AFGD{|{V)r)@?s?QWt68X?dGMjsh*ZkRU{HFR&chb z94Lq<$j4C({AY2DqBP!VN7jdvG>b6!7~c=W6eIk!egf5gU%ItVijEK z--Fa6f>j`Nu^_GZn>bX)m zYZH4yzIJt0I6L)T8UsnUiqb6-$&dZS{h2`Kt95lEur(C|lY{~xK1eY-i9yxJ|CPauu56CRe6MA8mg$j~dJS^FrGx{d1kO@|0aB2&l)npou6nyUwpg}L&01CK zd_IA;izq@WN3;zXWsDsS*mO!DCnV!{Smc7|r-M(8G=`VRCH5g;m027LX% z39?tbrDdJ({U^A=d)$Q2a3eMt%Tv_-awo1YMT_;Jx27uCqY6KwN<_5l)sz>44B<~Y zHlyabX7YMYBckULq$T1;iMW?DW4BCPwS - interaction.spring + model.spring.create diff --git a/assets/icons/model.spring.create.vectorcraft b/assets/icons/model.spring.create.vectorcraft new file mode 100644 index 0000000000000000000000000000000000000000..82a2de2c3c6cb76cd6fb9854eea49e374df4bc69 GIT binary patch literal 1896 zcmV-u2bcICiwFP!00000|K(b3Z`(Ey{wsz)3x-6!2iS)$E4rfHS|me@p$H6GqHQj+ zsOcq%7x=&LNQzFRC@w`O!+;uaVv!W@(cN>$<5Bn4id9?OSZCG^NgGzB6<#%#4Odl7 zSz(>|wzVb&sW5cHs4rsFIUhB1=E>ZFu^CI-Eh%7<`-l5n+iGZ&6S!x~GLp~hvZBR$ zo>l}m1nz1J+SJhdnck7ShK7{otXlBF(0Z|6LKCiH20a z46Fba7k}qiXU^QQt)2RD=UQjT7fs!-y=~PScDqOcG68|Oy0*?%IIoH7s90SuHkej* z*Rj8{tsL)2H3sA(?+_@iX_0Xh&93Af;tus-y!~@UY63VpgobR(oNGdf=IUwil54D~ zT4+c&G|wth97G~eNLdnG;UXn`J_?gaL{hS~qInJiOH8FTIv_R6;r>TJghimE5?Zdf zs^XA;BCxi&ibE9nfo(1MRpLgWA35QEzqb#v>ZW4XWJFrQ3Sz+%-qNhufTBR-SNkW* zDKT`~@*lL;?UMhLPXo`lVV$l^`0LJ*?V`B@R-nbhb&~k5ivm9hlfXwI^s=^5-5}Np z;us#n@f<&ht`q|6a?2P9<&VNA_kbI&Ge>H2G`P|ARkxho8_11S_xSYwJnmr^*CXKM zmH;X1a`_2O6dy>85&VtauTV%tt?xM!)`iq4Btiz0;J>oD=pZswGJ8g*z$VlAZg<4g-Fd;*#eJtS$KLe=%;!?($S` z65*AzQSm1nHbj6RdUOqH4Q#j9;thv{Plm+9X(dlZ%Ok>csR&pNzR+IZ;wIe?iMLBe zzWojcxjGVS);U|^{2}2fS)*pQ53_vT(8kP$95^?g_`j4fO;_sAIZm$)nd5YnaXPZ- zOT8eOq09vFU=gX6z?J+&9Rk-E?2Wa7d;;fNli3M)AQAxtzr*7N9N&p>OreBoSrkwL zxhg6shec8o%A4&{B+^DmV&8K;Wuq{0+$3~6Ac*&%AkYT#Cxm8pYm6jDUZ7%I;)##1 zbkcEYB9z&?iog!RFp=reC62hJjPZP=TT`;YYRj#nG!ym6g5Wo-vPc`qDH=gYI*mdi zJzD7mN%CkTdAo&ebDo_T(SnGlBaX^2kyM91Xa68LKM{iO^tgfCfV;PTD0%aDiwz;T z2R9kBNEh(5ft>Gfagxg@jf5h)sW=h}g?N)-nCc@Gwm=F;$A%6_w;9nl%s*{KP9=wk zuSof8=BY9tBbj?h$2s?OzBC7NlFH)BAI&j)( zzEGBv#t@A{;r2qQ@Y+aUA)ni-%#DM(5-ovj{6&iB4lG4L_!GLATS!_ zh#p56=u`EA)u!5rWu_56#}1jQ>05>trit<}Y3)k-_uh)#nufASy6;13Caygd zQ}q`(<6x@H5r{~E4vUIR)bztSTdjyGGX^=k$`K_krsnMD8wRp)W4!-Xidf53qqpkI zG)*7OB`!>JpQsl+`egtgP2}h~gUW{K_95(<$L|&UJedgTN3n~xHjU_m9NClcLlH&g zt=dSw;=E{$SC^o-lRtqtSg*R ze$-_@{5(6{Lp83+yX$cL{w)5<1LXe@C`n_bturCG@Y9H0Y-tmbx) z&plQgIX*Wpo5ALvA4SQ*dh7;$;CyzH2IU+-ss=CHq8`opUW;~`<&l@uqh=oT8p{2Y zlF>xr^>nMZS9>|-YVL<#PPtdA2EABUtSx(m=yw&~@l*Y=BK3!#5@$-Z2fOT`{aO%K zJCuL;DK?f8PmmHBONk~(35}(M6Ql&jQi2Imd}As81Sy`e6mNnQWGn?ukm4Fk`CjUC z)!WU{^`?2wYh|zT`Gk8|L|s?bhz<#($gzV7n@$1bgkt;y7O~(t@#IsZ4#P{-68jK* zSy`M4$pf9l>+n24NVtyv&;MJ{e`S>n{$ha z5wlA^!g{^;L`DxnU+SEE^t5yE7D@M#8h_1hv1E7M`v!y4*4d}F$sxtD{vYT6y3d7< zJlC%T-|S}ZMS<+~UINmu_x=~U@gM2NJE9xD!PuN)?w31peJQnQA9{DHLXTACKT#FW zXm{Hs%egXmJZWi0)bEOl7fEWKJr^J?Lfja2yUgef>-mass`, `-inertia`, `-spring`, made +//! unique against the model's sets): `MASS` with `*MASS` (one mass for every direction), +//! `ROTARYI` with `*ROTARY INERTIA` (I11, I22, I33, I12, I13, I23), and `SPRING1`/`DASHPOT1` to +//! the ground and `SPRING2`/`DASHPOT2` between two nodes along a global degree of freedom, or +//! `SPRINGA`/`DASHPOTA` along the line of the nodes, with `*SPRING`/`*DASHPOT`. A mass that +//! differs by direction and a spring along a skew fixed vector are reported (they need a local +//! orientation, not written yet). //! //! | ModelCraft | Plane stress | Plane strain | Axisymmetric | //! |---|---|---|---| @@ -137,13 +149,16 @@ pub fn write_resolved(model: &Model, r: &Resolved) -> Result { crate::refuse_fluid(FORMAT, model)?; let by_type = elements_by_type(model, r)?; check_expressible(model)?; + let discrete = discrete_groups(model, r)?; let mut d = Deck::default(); d.heading(model); d.nodes(model)?; d.elements(&by_type); + d.discrete_elements(&discrete); d.sets(model, r)?; d.materials(model)?; d.sections(model)?; + d.discrete_properties(&discrete)?; d.couplings(model, r); d.steps(model, r)?; Ok(d.out) @@ -198,6 +213,7 @@ pub fn element_type(kind: ElementKind, section: &Section) -> std::result::Result (S::Shell { .. }, K::Tri6, ..) => return Err("tri6 shells (Abaqus has no general-purpose six-node triangular shell)".into()), (S::Beam { .. }, K::Beam2, ..) => "B31", (S::Rod { .. }, K::Rod2, ..) => "T3D2", + (S::Rod { .. }, K::Rod3, ..) => "T3D3", (s, ..) => return Err(format!("{} elements in a {} section", kind.name(), s.name())), }; Ok(if f.hybrid { format!("{plain}H") } else { plain.to_string() }) @@ -216,7 +232,9 @@ fn elements_by_type<'m>(model: &'m Model, r: &Resolved) -> Result { count += elements.len(); by_type.entry(t).or_default().extend(elements); @@ -284,6 +302,90 @@ fn check_expressible(model: &Model) -> Result<()> { unsupported(FORMAT, bad) } +/// The property of a group of discrete elements. +enum DiscreteProperty { + Mass(f64), + Inertia([f64; 6]), + /// `*SPRING` or `*DASHPOT`: the degrees of freedom line (empty for the axial form) and the + /// stiffness or damping coefficient. + Spring { + dashpot: bool, + dofs: String, + value: f64, + }, +} + +/// The discrete elements of one definition: element type, element set, members and property. +struct DiscreteGroup { + element_type: &'static str, + elset: String, + members: Vec<(u32, Vec)>, + property: DiscreteProperty, +} + +/// Point masses, springs and dashpots as element groups (one per definition and element type, +/// labels counting up after the mesh's), or everything Abaqus cannot take. +fn discrete_groups<'r>(model: &Model, r: &'r Resolved) -> Result> { + let mut next = model.elements().keys().next_back().map_or(1, |l| l.saturating_add(1)); + let mut label = || -> Result { + let l = next; + if l > modelcraft_model::MAX_LABEL { + return Err(Error::Limit(format!( + "no element label left below {} for point masses, springs and dashpots", + modelcraft_model::MAX_LABEL as u64 + 1 + ))); + } + next += 1; + Ok(l) + }; + let mut taken: BTreeSet = model.sets().keys().map(|n| n.to_ascii_uppercase()).collect(); + let mut groups: Vec = Vec::new(); + // The (definition, element type) the last group holds. + let mut current: Option<(&str, &'static str)> = None; + let mut push = |name: &'r str, suffix: &str, element_type: &'static str, property: DiscreteProperty, member: (u32, Vec)| { + let key = (name, element_type); + if current.as_ref() == Some(&key) + && let Some(g) = groups.last_mut() + { + g.members.push(member); + return; + } + let elset = crate::unique_name(&format!("{name}-{suffix}"), &mut taken); + groups.push(DiscreteGroup { element_type, elset, members: vec![member], property }); + current = Some(key); + }; + let mut bad = BTreeSet::new(); + for m in &r.masses { + if m.mass.iter().any(|c| *c != 0.0) { + if m.mass[0] != m.mass[1] || m.mass[0] != m.mass[2] { + bad.insert(format!("point mass `{}`: a mass that differs by direction (needs a local orientation, not written yet)", m.name)); + } + push(&m.name, "mass", "MASS", DiscreteProperty::Mass(m.mass[0]), (label()?, vec![m.node])); + } + } + for m in &r.masses { + if m.inertia.iter().any(|c| *c != 0.0) { + push(&m.name, "inertia", "ROTARYI", DiscreteProperty::Inertia(m.inertia), (label()?, vec![m.node])); + } + } + for s in &r.springs { + let what = s.kind(); + let line = matches!(s.direction, modelcraft_model::ResolvedDirection::Line(_)); + let (element_type, dofs, nodes) = match (s.second, crate::axis_dof(s.direction), line) { + (Some(b), _, true) => (if s.dashpot { "DASHPOTA" } else { "SPRINGA" }, String::new(), vec![s.first, b]), + (Some(b), Some(d), false) => (if s.dashpot { "DASHPOT2" } else { "SPRING2" }, format!("{d}, {d}"), vec![s.first, b]), + (None, Some(d), _) => (if s.dashpot { "DASHPOT1" } else { "SPRING1" }, d.to_string(), vec![s.first]), + (_, None, _) => { + bad.insert(format!("{what} `{}`: a direction along no global axis (needs a local orientation, not written yet)", s.name)); + continue; + } + }; + push(&s.name, what, element_type, DiscreteProperty::Spring { dashpot: s.dashpot, dofs, value: s.value }, (label()?, nodes)); + } + unsupported(FORMAT, bad)?; + Ok(groups) +} + /// Append a real number in the shortest form that reads back to the same value: `0.1`, /// `210000.0`, `7.85E-9`, `1.0E300`. fn push_number(out: &mut String, v: f64) -> Result<()> { @@ -342,15 +444,7 @@ fn coupling_surfaces(model: &Model, r: &Resolved) -> BTreeMap { r.sets.iter().filter(|(n, s)| !s.faces.is_empty() || !r.planar_edges(n).is_empty()).map(|(n, _)| n.to_ascii_uppercase()).collect(); let mut out = BTreeMap::new(); for name in model.couplings().keys() { - let base = format!("{name}-nodes"); - let mut candidate = base.clone(); - let mut k = 2u32; - while taken.contains(&candidate.to_ascii_uppercase()) { - candidate = format!("{base}-{k}"); - k = k.saturating_add(1); - } - taken.insert(candidate.to_ascii_uppercase()); - out.insert(name.clone(), candidate); + out.insert(name.clone(), crate::unique_name(&format!("{name}-nodes"), &mut taken)); } out } @@ -581,6 +675,41 @@ impl Deck { } } + fn discrete_elements(&mut self, groups: &[DiscreteGroup]) { + if groups.is_empty() { + return; + } + self.comment("Point masses, springs and dashpots"); + for g in groups { + self.line(&format!("*ELEMENT, TYPE={}, ELSET={}", g.element_type, g.elset)); + for (label, nodes) in &g.members { + let nodes: Vec = nodes.iter().map(u32::to_string).collect(); + self.line(&format!("{label}, {}", nodes.join(", "))); + } + } + } + + fn discrete_properties(&mut self, groups: &[DiscreteGroup]) -> Result<()> { + for g in groups { + match &g.property { + DiscreteProperty::Mass(m) => { + self.line(&format!("*MASS, ELSET={}", g.elset)); + self.numbers_line("", &[*m])?; + } + DiscreteProperty::Inertia(i) => { + self.line(&format!("*ROTARY INERTIA, ELSET={}", g.elset)); + self.numbers_line("", i)?; + } + DiscreteProperty::Spring { dashpot, dofs, value } => { + self.line(&format!("{}, ELSET={}", if *dashpot { "*DASHPOT" } else { "*SPRING" }, g.elset)); + self.line(if dofs.is_empty() { "," } else { dofs }); + self.numbers_line("", &[*value])?; + } + } + } + Ok(()) + } + fn sets(&mut self, model: &Model, r: &Resolved) -> Result<()> { if r.sets.is_empty() { return Ok(()); diff --git a/crates/decks/src/lib.rs b/crates/decks/src/lib.rs index 076d236..973aaa9 100644 --- a/crates/decks/src/lib.rs +++ b/crates/decks/src/lib.rs @@ -16,7 +16,7 @@ pub mod nastran; use std::collections::{BTreeMap, BTreeSet}; use modelcraft_model::nodal::edge_pressure; -use modelcraft_model::{Diagnostic, Load, Model, Resolved}; +use modelcraft_model::{Diagnostic, Load, Model, Resolved, ResolvedDirection}; #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum Error { @@ -158,3 +158,28 @@ pub fn edge_pressure_loads<'a>(model: &Model, r: &Resolved, active: impl IntoIte out.retain(|_, v| *v != 0.0); Ok(out) } + +/// The global degree of freedom (1–6) a spring or dashpot direction acts along, when it is one +/// (a vector or line along an axis counts: the sign does not matter to a spring). +pub(crate) fn axis_dof(d: ResolvedDirection) -> Option { + match d { + ResolvedDirection::Dof(k) => Some(k), + ResolvedDirection::Vector(n) | ResolvedDirection::Line(n) => { + let axis = n.iter().position(|c| c.abs() == 1.0)?; + (n.iter().filter(|c| **c != 0.0).count() == 1).then_some(axis as u8 + 1) + } + } +} + +/// `base`, or `base-2`, `base-3`, … : the first that is not in `taken` (compared without case), +/// which it is then added to. +pub(crate) fn unique_name(base: &str, taken: &mut BTreeSet) -> String { + let mut candidate = base.to_string(); + let mut k = 2u32; + while taken.contains(&candidate.to_ascii_uppercase()) { + candidate = format!("{base}-{k}"); + k = k.saturating_add(1); + } + taken.insert(candidate.to_ascii_uppercase()); + candidate +} diff --git a/crates/decks/src/lsdyna/elements.rs b/crates/decks/src/lsdyna/elements.rs index 980dbc6..855edef 100644 --- a/crates/decks/src/lsdyna/elements.rs +++ b/crates/decks/src/lsdyna/elements.rs @@ -38,7 +38,7 @@ impl Keyword { ElementKind::Hex20 => Keyword::SolidH20, ElementKind::Tet4 | ElementKind::Tet10 | ElementKind::Wedge6 | ElementKind::Hex8 => Keyword::Solid, ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => Keyword::Shell, - ElementKind::Rod2 | ElementKind::Beam2 => Keyword::Beam, + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => Keyword::Beam, } } @@ -56,6 +56,7 @@ impl Keyword { pub(crate) fn order(kind: ElementKind) -> &'static [usize] { match kind { ElementKind::Rod2 | ElementKind::Beam2 => &[0, 1], + ElementKind::Rod3 => &[0, 1, 2], ElementKind::Tri3 => &[0, 1, 2, 2], ElementKind::Tri6 => &[0, 1, 2, 2, 3, 4, 5, 5], ElementKind::Quad4 => &[0, 1, 2, 3], @@ -141,6 +142,7 @@ pub(crate) fn elform(kind: ElementKind, section: &SectionKind, f: Formulation, s }, SectionKind::Rod { .. } => match kind { ElementKind::Rod2 => plain(3), + ElementKind::Rod3 => no("LS-DYNA's trusses have two nodes; there is no counterpart for"), _ => no("a rod section cannot hold"), }, SectionKind::Beam { .. } => no("beam sections are not written to LS-DYNA decks yet (planned with beam profiles):"), diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 8e6cd6a..c011f97 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -31,6 +31,15 @@ //! (with `*LOAD_BODY_PARTS` when it acts on a set); output as `*DATABASE_BINARY_D3PLOT`, //! `*DATABASE_SPCFORC` and `*DATABASE_EXTENT_BINARY`. //! +//! Point masses are `*ELEMENT_MASS` (one mass for every direction) or `*ELEMENT_INERTIA` (the +//! inertia tensor's components, with the mass); springs and dashpots are `*ELEMENT_DISCRETE` with +//! one `*PART`, `*SECTION_DISCRETE` (torsional for a rotational degree of freedom) and +//! `*MAT_SPRING_ELASTIC` or `*MAT_DAMPER_VISCOUS` per definition, acting along the line of their +//! nodes or along a `*DEFINE_SD_ORIENTATION` vector. A spring or dashpot to the ground connects +//! to an extra node at the same place, held by `*BOUNDARY_SPC_NODE` (labels after the mesh's). +//! Element, part, section and material ids continue after the mesh's. A mass that differs by +//! direction and the three-node truss (rod3) are reported. +//! //! Several steps are independent analyses in the model (nothing carries over), which is what //! LS-DYNA's `*CASE` runs: each step is a case between `*CASE_BEGIN_n` and `*CASE_END_n`, run in //! order when LS-DYNA is started with `CASE` on its command line. A single step is written @@ -54,7 +63,7 @@ use std::collections::{BTreeMap, BTreeSet}; use modelcraft_model::{ BoundaryCondition, CouplingKind, Diagnostic, Dofs, Elastic, Element, ElementKind, Hardening, Load, Material, Model, OutputVariable, Plastic, - Procedure, Resolved, Section, SectionKind, Step, Target, + Procedure, Resolved, ResolvedDirection, Section, SectionKind, Step, Target, }; use crate::{Error, Result}; @@ -161,6 +170,7 @@ fn write_with(model: &Model, r: &Resolved) -> Result { element_cards(&mut d, &plan)?; sets(&mut d, model, r, &plan)?; couplings(&mut d, model, r, &plan)?; + discrete(&mut d, model, r, &plan)?; if let Some(parts) = &plan.gravity_parts { d.line("*SET_PART_LIST"); d.header(&["sid"]); @@ -312,6 +322,16 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result> { None } }; + if r.is_planar() { + for name in model.sections().iter().filter(|(_, s)| matches!(s.kind, SectionKind::Rod { .. })).map(|(n, _)| n) { + refused.insert(format!("section `{name}`: trusses in a 2D model (LS-DYNA's 2D elements take no truss)")); + } + } + for m in &r.masses { + if m.mass[0] != m.mass[1] || m.mass[0] != m.mass[2] { + refused.insert(format!("point mass `{}`: a mass that differs by direction (LS-DYNA's mass elements take one mass)", m.name)); + } + } let mut curves = BTreeMap::new(); for name in materials.keys() { let Some(m) = model.material(name) else { continue }; @@ -646,6 +666,122 @@ fn couplings(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan) -> Result<( Ok(()) } +/// Point masses, springs and dashpots (see the module docs). +fn discrete(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan) -> Result<()> { + if r.masses.is_empty() && r.springs.is_empty() { + return Ok(()); + } + let limit = |what: &str| Error::Limit(format!("no {what} id left for point masses, springs and dashpots")); + let mut eid = model.elements().keys().next_back().copied().unwrap_or(0); + let mut next_eid = || -> Result { + eid = eid.checked_add(1).filter(|e| *e <= modelcraft_model::MAX_LABEL).ok_or_else(|| limit("element"))?; + Ok(eid) + }; + d.comment("Point masses"); + for m in &r.masses { + let mass = m.mass[0]; + if m.inertia.iter().all(|c| *c == 0.0) { + d.line("*ELEMENT_MASS"); + d.header(&["eid", "nid", "mass"]); + d.card([Field::from(next_eid()?), Field::from(m.node), Field::from(mass)]); + } else { + let [i11, i22, i33, i12, i13, i23] = m.inertia; + d.line("*ELEMENT_INERTIA"); + d.header(&["eid", "nid", "csid"]); + d.card([Field::from(next_eid()?), Field::from(m.node)]); + d.header(&["ixx", "ixy", "ixz", "iyy", "iyz", "izz", "mass"]); + d.card([i11, i12, i13, i22, i23, i33, mass].map(Field::from)); + } + } + if r.springs.is_empty() { + return Ok(()); + } + // Ids after the mesh's parts, the rigid bodies' parts and the materials. + let kinematic = model.couplings().values().filter(|c| c.kind == CouplingKind::Kinematic).count(); + let mut pid = u32::try_from(plan.parts.len() + kinematic).map_err(|_| limit("part"))?; + let mut mid = u32::try_from(plan.materials.len()).map_err(|_| limit("material"))?; + let mut vid = 0u32; + // Ground nodes: labels after the mesh's, skipping the reference points'. + let taken: BTreeSet = model.reference_points().values().map(|p| p.label).collect(); + let mut node = model.nodes().keys().next_back().copied().unwrap_or(0); + let mut ground: Vec<(u32, [f64; 3])> = Vec::new(); + let mut elements: Vec> = Vec::new(); + let mut current: Option<&str> = None; + let mut vid_of_part = 0u32; + for s in &r.springs { + let rotational = s.direction.is_rotational(); + if current != Some(s.name.as_str()) { + current = Some(&s.name); + pid = pid.checked_add(1).ok_or_else(|| limit("part"))?; + mid = mid.checked_add(1).ok_or_else(|| limit("material"))?; + d.comment(&format!("{} `{}`", s.kind(), s.name)); + if s.dashpot { + d.line("*MAT_DAMPER_VISCOUS_TITLE"); + d.line(&s.name); + d.header(&["mid", "dc"]); + } else { + d.line("*MAT_SPRING_ELASTIC_TITLE"); + d.line(&s.name); + d.header(&["mid", "k"]); + } + d.card([Field::from(mid), Field::from(s.value)]); + d.line("*SECTION_DISCRETE"); + d.header(&["secid", "dro"]); + d.card([Field::from(pid), Field::from(i32::from(rotational))]); + d.line("*PART"); + d.line(&s.name); + d.header(&["pid", "secid", "mid"]); + d.card([Field::from(pid), Field::from(pid), Field::from(mid)]); + let axis = (!matches!(s.direction, ResolvedDirection::Line(_))).then(|| s.direction.axis()); + let v = match axis { + Some([x, y, z]) => { + vid += 1; + d.line("*DEFINE_SD_ORIENTATION"); + d.header(&["vid", "iop", "xt", "yt", "zt"]); + d.card([Field::from(vid), Field::from(0), Field::from(x), Field::from(y), Field::from(z)]); + vid + } + None => 0, + }; + vid_of_part = v; + } + let second = match s.second { + Some(b) => b, + None => { + node = node.checked_add(1).ok_or_else(|| limit("node"))?; + while taken.contains(&node) { + node = node.checked_add(1).ok_or_else(|| limit("node"))?; + } + if node > modelcraft_model::MAX_LABEL { + return Err(limit("node")); + } + ground.push((node, model.label_xyz(s.first).ok_or_else(|| broken(format!("node {}", s.first)))?)); + node + } + }; + elements.push(vec![Field::from(next_eid()?), Field::from(pid), Field::from(s.first), Field::from(second), Field::from(vid_of_part)]); + } + d.line("*ELEMENT_DISCRETE"); + d.header(&["eid", "pid", "n1", "n2", "vid"]); + for e in elements { + d.card(e); + } + if !ground.is_empty() { + d.comment("The ground ends of springs and dashpots, held in every direction"); + d.line("*NODE"); + d.header(&["nid", "x", "y", "z"]); + for (l, [x, y, z]) in &ground { + d.card([Field::from(*l), Field::from(*x), Field::from(*y), Field::from(*z)]); + } + d.line("*BOUNDARY_SPC_NODE"); + d.header(&["nid", "cid", "dofx", "dofy", "dofz", "dofrx", "dofry", "dofrz"]); + for (l, _) in &ground { + d.card([Field::from(*l), Field::Blank].into_iter().chain([1; 6].map(Field::from))); + } + } + Ok(()) +} + /// LS-DYNA's degree-of-freedom code for prescribed motion and nodal loads: 1–3 translations, /// 5–7 rotations. fn dof_code(dof: u8) -> i32 { diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs index 297f2a7..3455abd 100644 --- a/crates/decks/src/nastran/mod.rs +++ b/crates/decks/src/nastran/mod.rs @@ -35,7 +35,13 @@ //! (`PPLANE`/`PLPLANE` with `CPLSTN*`/`CQUADX`, `CTRIAX6`, `PAXI`) differ between Nastran //! dialects, so a deck could not be read the same everywhere; //! - centrifugal loads: `RFORCE` about an axis through the origin of the basic system (speed in -//! revolutions per time unit). +//! revolutions per time unit); +//! - point masses: `CONM1`, the node's 6 × 6 mass matrix (the masses per direction on the +//! translations, the inertia tensor's components on the rotations, so no product-of-inertia sign +//! convention is involved); springs and dashpots along a global degree of freedom: `CELAS2` and +//! `CDAMP2` (the second grid left blank for one to the ground). Their ids follow the rigid +//! elements'. A spring or dashpot along a skew direction and the three-node truss (rod3) have +//! no counterpart here and are reported. //! //! Field widths, continuations and number forms are in [`card`]. Model content Nastran cannot //! express is reported as [`Error::Unsupported`], every item at once, and nothing is written. @@ -96,7 +102,7 @@ pub fn order(kind: ElementKind) -> &'static [usize] { /// The connection entry's name for an element kind. pub fn element_card(kind: ElementKind) -> &'static str { match kind { - ElementKind::Rod2 => "CROD", + ElementKind::Rod2 | ElementKind::Rod3 => "CROD", ElementKind::Beam2 => "CBEAM", ElementKind::Tri3 => "CTRIA3", ElementKind::Tri6 => "CTRIA6", @@ -257,6 +263,8 @@ struct Deck<'m> { element_pid: BTreeMap, /// Rigid elements: (id, coupling name, coupling). rigid: Vec<(u32, &'m str, &'m Coupling)>, + /// Point masses, springs and dashpots: (comment, card). + discrete: Vec<(String, Card)>, /// Node → (dependent degrees of freedom, coupling that makes them dependent). dependent: BTreeMap, /// Constraint sets by id. @@ -304,6 +312,7 @@ impl<'m> Deck<'m> { properties: Vec::new(), element_pid: BTreeMap::new(), rigid: Vec::new(), + discrete: Vec::new(), dependent: BTreeMap::new(), spc_sets: Vec::new(), load_sids: BTreeMap::new(), @@ -315,6 +324,7 @@ impl<'m> Deck<'m> { self.plan_materials()?; self.plan_sections()?; self.plan_couplings()?; + self.plan_discrete()?; self.plan_constraints()?; self.plan_loads()?; for sub in &mut self.subcases { @@ -447,6 +457,9 @@ impl<'m> Deck<'m> { self.add_property(sname, Property::Beam { mid, profile, orientation }, elements.iter().copied())?; } SectionKind::Rod { area } => { + if elements.iter().any(|e| kind_of(*e) == Some(ElementKind::Rod3)) { + self.problems.insert(format!("section `{sname}`: three-node trusses (rod3; Nastran's rods have two grids)")); + } let Some(mid) = ids.isotropic else { continue }; self.add_property(sname, Property::Rod { mid, area }, elements.iter().copied())?; } @@ -520,6 +533,44 @@ impl<'m> Deck<'m> { Ok(()) } + /// Point masses, springs and dashpots, numbered after the rigid elements. + fn plan_discrete(&mut self) -> Result<()> { + let r = self.r; + let mut next = self.rigid.last().map(|e| e.0).or_else(|| self.model.elements().keys().next_back().copied()).unwrap_or(0); + let mut id = || -> Result { + next = next + .checked_add(1) + .filter(|n| *n <= MAX_LABEL) + .ok_or_else(|| Error::Limit(format!("no element id left below {} for point masses, springs and dashpots", MAX_LABEL as u64 + 1)))?; + Ok(next) + }; + for m in &r.masses { + let [mx, my, mz] = m.mass; + let [i11, i22, i33, i12, i13, i23] = m.inertia; + // M11, M21, M22, M31, M32, M33, then the rotations' rows: M41…M44, M51…M55, M61…M66. + let card = Card::new("CONM1") + .int(id()?) + .int(m.node) + .blank() + .reals([mx, 0.0, my, 0.0, 0.0, mz, 0.0, 0.0, 0.0, i11, 0.0, 0.0, 0.0, i12, i22, 0.0, 0.0, 0.0, i13, i23, i33]); + self.discrete.push((format!("point mass {}", m.name), card)); + } + for s in &r.springs { + let what = s.kind(); + let Some(dof) = crate::axis_dof(s.direction) else { + self.problems.insert(format!("{what} `{}`: a direction along no global axis (CELAS and CDAMP act along grid components)", s.name)); + continue; + }; + let card = Card::new(if s.dashpot { "CDAMP2" } else { "CELAS2" }).int(id()?).real(s.value).int(s.first).int(dof); + let card = match s.second { + Some(b) => card.int(b).int(dof), + None => card, + }; + self.discrete.push((format!("{what} {}", s.name), card)); + } + Ok(()) + } + /// One constraint set per distinct step constraint, checked against coupled nodes. fn plan_constraints(&mut self) -> Result<()> { let (model, r) = (self.model, self.r); @@ -780,6 +831,13 @@ impl<'m> Deck<'m> { }; write_named(out, &format!("coupling {name}"), *eid, card)?; } + if !self.discrete.is_empty() { + comment(out, "point masses, springs and dashpots"); + } + for (what, card) in &self.discrete { + comment(out, what); + card.write(out)?; + } Ok(()) } diff --git a/crates/decks/tests/discrete.rs b/crates/decks/tests/discrete.rs new file mode 100644 index 0000000..1c2aa3e --- /dev/null +++ b/crates/decks/tests/discrete.rs @@ -0,0 +1,166 @@ +//! Trusses (rod3), point masses, rotary inertia, springs and dashpots in every deck format: +//! what each writes, and what each reports as unsupported. + +use modelcraft_decks::{Error, abaqus, lsdyna, nastran}; +use modelcraft_model::{ + BoundaryCondition, Direction, Elastic, ElementKind, Formulation, Mass, Material, Model, PointMass, Procedure, ReferencePoint, Section, + SectionKind, Set, Spring, Step, Target, UnitSystem, +}; + +/// A rod2 (or rod3) bar from node 1 to node 2 and a rod2 on to node 4, a reference point beyond +/// it; a motor mass at the tip, a hub with inertia on the point, a spring to the ground across +/// the tip, a spring along the line from the tip to the hub, a dashpot to the ground. +fn model(rod3: bool) -> Model { + let mut m = Model::default(); + m.set_units(UnitSystem::MmTonneS); + for (l, x) in [(1, 0.0), (2, 1000.0), (3, 500.0), (4, 2000.0)] { + m.set_node(l, [x, 0.0, 0.0]).unwrap(); + } + if rod3 { + m.set_element(1, ElementKind::Rod3, vec![1, 2, 3]).unwrap(); + } else { + m.set_element(1, ElementKind::Rod2, vec![1, 2]).unwrap(); + } + m.set_element(2, ElementKind::Rod2, vec![2, 4]).unwrap(); + m.add_set("bars", Set { elements: vec![1, 2], ..Set::default() }).unwrap(); + m.add_set("base", Set { nodes: vec![1], ..Set::default() }).unwrap(); + m.add_set("tip", Set { nodes: vec![4], ..Set::default() }).unwrap(); + m.add_material( + "Steel", + Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() }, + ) + .unwrap(); + m.add_section( + "bars", + Section { set: "bars".into(), material: "Steel".into(), kind: SectionKind::Rod { area: 20.0 }, formulation: Formulation::default() }, + ) + .unwrap(); + m.add_reference_point("Hub", ReferencePoint { xyz: [2500.0, 0.0, 0.0], label: 0 }).unwrap(); + let tip = Target::Set("tip".into()); + let hub = Target::ReferencePoint("Hub".into()); + m.add_mass("motor", PointMass { target: tip.clone(), mass: Some(Mass::Isotropic(0.002)), inertia: None }).unwrap(); + m.add_mass("hub", PointMass { target: hub.clone(), mass: Some(Mass::Isotropic(0.01)), inertia: Some([1.0, 2.0, 3.0, 0.1, 0.0, 0.0]) }).unwrap(); + m.add_spring("mount", Spring::Spring { first: tip.clone(), second: None, direction: Direction::Dof(2), stiffness: 50.0 }).unwrap(); + m.add_spring("link", Spring::Spring { first: tip.clone(), second: Some(hub), direction: Direction::Line, stiffness: 1000.0 }).unwrap(); + m.add_spring("damper", Spring::Dashpot { first: tip, second: None, direction: Direction::Vector([0.0, 0.0, -2.0]), coefficient: 0.5 }).unwrap(); + m.add_boundary_condition("fix", BoundaryCondition::Fixed { target: Target::Set("base".into()) }).unwrap(); + let mut step = Step::new("Modes", Procedure::Frequency { modes: 3, min_hz: None, max_hz: None }); + step.boundary_conditions.insert("fix".into()); + m.add_step(step, None).unwrap(); + m +} + +fn unsupported(e: std::result::Result) -> String { + match e { + Err(Error::Unsupported { what, .. }) => what, + other => panic!("expected unsupported, got {other:?}"), + } +} + +#[test] +fn abaqus_writes_every_discrete_element_and_rod3() { + let deck = abaqus::write(&model(true)).unwrap(); + for expected in [ + "*ELEMENT, TYPE=T3D3\n1, 1, 2, 3\n", + "*ELEMENT, TYPE=MASS, ELSET=hub-mass\n3, 99999999\n", + "*ELEMENT, TYPE=MASS, ELSET=motor-mass\n4, 4\n", + "*ELEMENT, TYPE=ROTARYI, ELSET=hub-inertia\n5, 99999999\n", + "*ELEMENT, TYPE=DASHPOT1, ELSET=damper-dashpot\n6, 4\n", + "*ELEMENT, TYPE=SPRINGA, ELSET=link-spring\n7, 4, 99999999\n", + "*ELEMENT, TYPE=SPRING1, ELSET=mount-spring\n8, 4\n", + "*MASS, ELSET=hub-mass\n0.01\n", + "*ROTARY INERTIA, ELSET=hub-inertia\n1.0, 2.0, 3.0, 0.1, 0.0, 0.0\n", + "*DASHPOT, ELSET=damper-dashpot\n3\n0.5\n", + "*SPRING, ELSET=link-spring\n,\n1000.0\n", + "*SPRING, ELSET=mount-spring\n2\n50.0\n", + ] { + assert!(deck.contains(expected), "{expected}\n{deck}"); + } + // An element set name that a model set already has is made unique. + let mut m = model(false); + m.add_set("motor-mass", Set { nodes: vec![4], ..Set::default() }).unwrap(); + assert!(abaqus::write(&m).unwrap().contains("ELSET=motor-mass-2\n")); + // Skew directions and masses per direction need orientations. + let mut m = model(false); + m.replace_spring( + "mount", + Spring::Spring { first: Target::Set("tip".into()), second: None, direction: Direction::Vector([0.0, 1.0, 1.0]), stiffness: 1.0 }, + ) + .unwrap(); + m.replace_mass("motor", PointMass { target: Target::Set("tip".into()), mass: Some(Mass::Directional([1.0, 2.0, 3.0])), inertia: None }).unwrap(); + let what = unsupported(abaqus::write(&m)); + assert!(what.contains("spring `mount`") && what.contains("point mass `motor`"), "{what}"); +} + +#[test] +fn nastran_writes_conm1_celas2_cdamp2() { + let deck = nastran::write(&model(false)).unwrap(); + eprintln!("{deck}"); + for expected in ["CONM1*", "CELAS2*", "CDAMP2*", "$ point mass hub", "$ spring link", "$ dashpot damper"] { + assert!(deck.contains(expected), "{expected}\n{deck}"); + } + // The hub's 6 × 6 mass matrix: mass on the translations, the inertia tensor on the rotations. + let conm1: Vec<&str> = deck.lines().skip_while(|l| !l.starts_with("CONM1*")).take(6).collect(); + let fields: Vec = conm1 + .iter() + .flat_map(|l| l.get(8..).unwrap_or("").as_bytes().chunks(16).map(|c| std::str::from_utf8(c).unwrap().trim().to_string())) + .filter(|t| !t.is_empty()) + .filter_map(|t| t.parse().ok()) + .collect(); + assert_eq!(fields.first().copied(), Some(3.0), "element id after the mesh's: {conm1:?}"); + assert!( + fields.ends_with(&[0.01, 0.0, 0.01, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.1, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0]), + "{fields:?}" + ); + // The line spring runs along X: component 1 at both grids; the dashpot along −Z: component 3. + let link = deck.lines().skip_while(|l| !l.contains("spring link")).nth(1).unwrap().to_string(); + assert!(link.starts_with("CELAS2*") && link.contains("1000.") && link.ends_with("4 1"), "{link}"); + let rest = deck.lines().skip_while(|l| !l.contains("spring link")).nth(2).unwrap().to_string(); + assert!(rest.contains("99999999") && rest.trim_end().ends_with('1'), "{rest}"); + let what = unsupported(nastran::write(&model(true))); + assert!(what.contains("rod3"), "{what}"); + let mut m = model(false); + m.replace_spring( + "link", + Spring::Spring { first: Target::Set("tip".into()), second: None, direction: Direction::Vector([1.0, 1.0, 0.0]), stiffness: 1.0 }, + ) + .unwrap(); + assert!(unsupported(nastran::write(&m)).contains("spring `link`")); +} + +#[test] +fn lsdyna_writes_masses_inertia_and_discrete_springs() { + let deck = lsdyna::write(&model(false)).unwrap().text; + eprintln!("{deck}"); + for expected in [ + "*ELEMENT_MASS", + "*ELEMENT_INERTIA", + "*MAT_SPRING_ELASTIC_TITLE\nlink\n", + "*MAT_DAMPER_VISCOUS_TITLE\ndamper\n", + "*SECTION_DISCRETE", + "*DEFINE_SD_ORIENTATION", + "*ELEMENT_DISCRETE", + "*BOUNDARY_SPC_NODE", + ] { + assert!(deck.contains(expected), "{expected}\n{deck}"); + } + // Two grounded links (damper, mount) get two ground nodes after the mesh's (5, 6); the + // line spring has no orientation vector (vid 0, along its nodes). + let discrete: Vec> = deck + .lines() + .skip_while(|l| *l != "*ELEMENT_DISCRETE") + .skip(1) + .take_while(|l| !l.starts_with('*')) + .filter(|l| !l.starts_with('$')) + .map(|l| l.split_whitespace().map(|f| f.parse().unwrap()).collect()) + .collect(); + assert_eq!(discrete.len(), 3); + assert_eq!(discrete[0][2..4], [4, 5], "damper to its ground node"); + assert_eq!(discrete[1][2..4], [4, 99_999_999], "link to the hub"); + assert_eq!(discrete[1].get(4).copied().unwrap_or(0), 0, "along the line"); + assert_eq!(discrete[2][2..4], [4, 6], "mount to its ground node"); + assert!(unsupported(lsdyna::write(&model(true))).contains("rod3")); + let mut m = model(false); + m.replace_mass("motor", PointMass { target: Target::Set("tip".into()), mass: Some(Mass::Directional([1.0, 2.0, 3.0])), inertia: None }).unwrap(); + assert!(unsupported(lsdyna::write(&m)).contains("point mass `motor`")); +} diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index caac4a0..abef15d 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -306,8 +306,22 @@ pub static CATALOGUE: &[Capability] = &[ soon("interaction.contact.find", "Find Contacts", Interactions, "Contact", "M7", "Detect contact pairs automatically"), soon("interaction.rigid_wall", "Rigid Wall", Interactions, "Contact", "M10", "Rigid walls for explicit dynamics"), soon("interaction.bolt", "Bolt Load", Interactions, "Connectors", "M7", "Bolt pretension"), - soon("interaction.spring", "Spring", Interactions, "Connectors", "M7", "Springs and dashpots"), - soon("interaction.mass", "Point Mass", Interactions, "Connectors", "M7", "Point masses and inertia"), + ready_more( + "model.spring.create", + "Spring", + Interactions, + "Connectors", + "Spring or dashpot to the ground or between two points, along a direction or their line", + &["model.spring.edit", "model.spring.delete"], + ), + ready_more( + "model.mass.create", + "Point Mass", + Interactions, + "Connectors", + "Point mass and rotary inertia at nodes or a reference point", + &["model.mass.edit", "model.mass.delete"], + ), soon("interaction.connector", "Connector", Interactions, "Connectors", "M7", "Connector elements between points"), soon("interaction.fsi", "FSI Interface", Interactions, "Coupled", "M13", "The fluid–structure interface"), // ---- Solve ---- diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs index c95217a..a5c913b 100644 --- a/crates/engine/src/cmd/define.rs +++ b/crates/engine/src/cmd/define.rs @@ -1,13 +1,14 @@ //! Model definition commands: create, edit and delete sets, materials, sections, reference -//! points, couplings, boundary conditions, loads, fluid domains and steps, and the step activity -//! table. +//! points, couplings, point masses, springs and dashpots, boundary conditions, loads, fluid +//! domains and steps, and the step activity table. //! //! Each definition's parameters are its JSON form in the model (the same form a `.modelcraft` //! file stores), plus its `name`. `create` fails when the name is taken, `edit` replaces an //! existing definition, `delete` refuses while another definition refers to it. use modelcraft_model::{ - BoundaryCondition, Coupling, FaceCondition, FluidBoundary, FluidDomain, Load, Material, Model, ReferencePoint, Section, Set, Step, + BoundaryCondition, Coupling, FaceCondition, FluidBoundary, FluidDomain, Load, Material, Model, PointMass, ReferencePoint, Section, Set, Spring, + Step, }; use serde::de::DeserializeOwned; use serde_json::{Value, json}; @@ -176,6 +177,34 @@ definition_commands!( COUPLING_DELETE ); +definition_commands!( + "point mass", + PointMass, + "mass", + "Point Mass", + r#"{"name": "motor", "target": {"set": "…"} | {"reference_point": "…"}, "mass": 2.5 | [mx, my, mz] (per global direction), "inertia": [I11, I22, I33, I12, I13, I23] (tensor about the node in global axes; needs rotations: a reference point)} — at every node of the target, each with the full value"#, + add_mass, + replace_mass, + remove_mass, + MASS_CREATE, + MASS_EDIT, + MASS_DELETE +); + +definition_commands!( + "spring", + Spring, + "spring", + "Spring or Dashpot", + r#"{"name": "mount", "type": "spring", "stiffness": 1000 | "type": "dashpot", "coefficient": 5, "first": {"set": "…"} | {"reference_point": "…"}, "second": (optional; without it each node of `first` connects to the ground; with it each end is one node), "direction": {"dof": 1..6} (4–6: torsional) | {"vector": [x, y, z]} | "line" (along the two ends, two-point only)} — dashpots act in dynamic steps; static and frequency steps ignore them"#, + add_spring, + replace_spring, + remove_spring, + SPRING_CREATE, + SPRING_EDIT, + SPRING_DELETE +); + definition_commands!( "boundary condition", BoundaryCondition, diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index f23d4eb..bc01e75 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -54,6 +54,12 @@ pub static ALL: &[CommandSpec] = &[ define::COUPLING_CREATE, define::COUPLING_EDIT, define::COUPLING_DELETE, + define::MASS_CREATE, + define::MASS_EDIT, + define::MASS_DELETE, + define::SPRING_CREATE, + define::SPRING_EDIT, + define::SPRING_DELETE, define::BOUNDARY_CONDITION_CREATE, define::BOUNDARY_CONDITION_EDIT, define::BOUNDARY_CONDITION_DELETE, diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index 54ce8da..f27423d 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -950,4 +950,32 @@ pub(crate) mod tests { // Defining steps edited the model, so the impact result is stale. assert_eq!(s.solved_steps().count(), 0); } + + /// A mass on a spring and a dashpot, defined, solved and written through commands; editing + /// and deleting them goes through the same checks as every definition. + #[test] + fn spring_mass_through_commands() { + let mut s = Session::default(); + let script = json!([ + {"command": "mesh.node.add", "params": {"label": 1, "xyz": [0, 0, 0]}}, + {"command": "model.set.create", "params": {"name": "n", "nodes": [1]}}, + {"command": "model.mass.create", "params": {"name": "m", "target": {"set": "n"}, "mass": 4}}, + {"command": "model.spring.create", "params": {"name": "k", "type": "spring", "first": {"set": "n"}, "direction": {"dof": 1}, "stiffness": 400}}, + {"command": "model.spring.create", "params": {"name": "c", "type": "dashpot", "first": {"set": "n"}, "direction": {"vector": [1, 0, 0]}, "coefficient": 1}}, + {"command": "model.boundary_condition.create", "params": {"name": "yz", "type": "displacement", "target": {"set": "n"}, "values": [null, 0, 0, null, null, null]}}, + {"command": "model.step.create", "params": {"name": "Modes", "procedure": {"type": "frequency", "modes": 1}, "boundary_conditions": ["yz"]}}, + ]); + s.run_script(&script).unwrap(); + let r = s.execute("solve.frequency", &json!({})).unwrap(); + let f = r["modes"][0]["frequency"].as_f64().unwrap(); + assert!((f - 10.0 / (2.0 * std::f64::consts::PI)).abs() < 1e-12, "{r}"); + assert!(r["warnings"].as_array().unwrap().iter().any(|w| w.as_str().unwrap().contains("dashpots `c`")), "{r}"); + assert!(s.execute("model.mass.edit", &json!({"name": "m", "target": {"set": "n"}, "mass": -1})).is_err()); + assert!(s.execute("model.set.delete", &json!({"name": "n"})).is_err(), "the set is in use"); + s.execute("model.mass.edit", &json!({"name": "m", "target": {"set": "n"}, "mass": 1})).unwrap(); + let r = s.execute("solve.frequency", &json!({})).unwrap(); + assert!((r["modes"][0]["frequency"].as_f64().unwrap() - 20.0 / (2.0 * std::f64::consts::PI)).abs() < 1e-12); + s.execute("model.spring.delete", &json!({"name": "c"})).unwrap(); + assert!(s.model.springs().len() == 1 && s.model.mass("m").is_some()); + } } diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index c25fc74..9773773 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -404,6 +404,13 @@ mod tests { json!({"name": "Bend", "procedure": {"type": "static"}, "loads": "push"}), json!({"step": "Bend", "loads": ["nope"], "active": "yes"}), json!({"step": "Bend", "boundary_conditions": [], "loads": []}), + json!({"name": "A", "target": {"set": "end"}, "mass": [1e308, -1, f64::MAX], "inertia": [1, 1, 1, 1e300, 0, 0]}), + json!({"name": "A", "target": {"reference_point": "RP"}, "inertia": [1, 2]}), + json!({"name": "A", "target": {"set": "end"}, "mass": "heavy"}), + json!({"name": "A", "type": "spring", "first": {"set": "end"}, "second": {"set": "end"}, "direction": "line", "stiffness": 1e308}), + json!({"name": "A", "type": "spring", "first": {"set": "end"}, "direction": {"dof": 255}, "stiffness": 1}), + json!({"name": "A", "type": "dashpot", "first": {"set": "end"}, "direction": {"vector": [1e-320, 0, 0]}, "coefficient": -1}), + json!({"name": "A", "type": "dashpot", "first": {"reference_point": "RP"}, "direction": {"dof": 4}, "coefficient": 1e300}), json!({"system": "furlongs"}), // Fluid domains, boundaries and steps. json!({"name": "A", "material": "Steel", "min": [0, 0, 0], "max": [1e308, 1, 1], "spacing": 1e-300}), diff --git a/crates/explicit/src/discrete.rs b/crates/explicit/src/discrete.rs new file mode 100644 index 0000000..51034b1 --- /dev/null +++ b/crates/explicit/src/discrete.rs @@ -0,0 +1,316 @@ +//! Trusses, point masses, springs and dashpots in explicit dynamics. +//! +//! - **Trusses** (rod2, rod3), elastic, with lumped (HRZ) mass. With nonlinear geometry a rod2 +//! is corotational: its strain is (l − L)/L from its current length l and its force acts along +//! its current axis (Belytschko et al. 2014, §2.4 and §4.5); a rod3 then is not supported yet. +//! Without it the internal force is K u with the kernel's stiffness ([`modelcraft_fem::truss`]). +//! The stress is the uniaxial tensor along the axis (the mean of the integration points'), and +//! the critical increment comes from the kernel's frequency bound (exact for rod2: L/c). +//! - **Point masses** add to the lumped mass of their node, per direction. +//! - **Springs** add the force k (cᵀ Δu) c and **dashpots** c_d (cᵀ Δv) c, with Δu, Δv the +//! relative motion of the ends along the unit direction c (a fixed direction; one along the line +//! of the ends follows that line as it turns under nonlinear geometry, with the stretch +//! l − L₀); dashpots use the velocity of the +//! previous half increment, v^{n−½}, as the bulk viscosity does (Belytschko, Liu, Moran, +//! Elkhodary 2014, §6.2). Spring work counts as internal energy, dashpot work as viscous +//! (dissipated) energy. +//! - **Stable increment** of springs and dashpots, per degree of freedom: with +//! ω² ≤ Σⱼ |Kᵢⱼ| / √(mᵢ mⱼ) (Gershgorin's bound on M^{-1/2} K M^{-1/2}, exact for one spring on +//! one mass) and the same row sum D of the damping, the central-difference scheme with lagged +//! damping is stable for Δt ≤ (2/ω)(√(1 + ξ²) − ξ) with ξ = D / (2ω) (Belytschko et al. 2014, +//! §6.3), which is 2/D for a dashpot alone. Rotary inertia and rotational springs need +//! rotational degrees of freedom, which this solver does not have yet. + +use modelcraft_fem::DofKind; +use modelcraft_fem::dense::Matrix; +use modelcraft_fem::truss::{Truss, uniaxial}; +use modelcraft_model::{ElementKind, Resolved, ResolvedDirection}; + +use crate::system::NODE_DOFS; +use crate::{Error, Result}; + +/// A truss element. +#[derive(Debug, Clone)] +pub struct Rod { + pub label: u32, + /// Position in [`crate::system::System::elements`]. + pub index: usize, + pub nodes: Vec, + pub material: usize, + pub young: f64, + /// Stiffness on the element's translations, node by node. + pub stiffness: Matrix, + /// Per integration point: the axial strain operator and the tangent. + pub strain: Vec<(Vec, [f64; 3])>, + /// The lumped mass per node for unit density. + pub nodal_volume: Vec, + /// The stable length: the critical increment is `length / wave_speed`. + pub length: f64, + /// E A and the rest length (a rod2's corotational force). + pub axial: f64, + pub rest_length: f64, +} + +/// What a truss is made of: Young's modulus, density and cross-sectional area. +#[derive(Debug, Clone, Copy)] +pub struct RodSection { + pub young: f64, + pub density: f64, + pub area: f64, +} + +impl Rod { + /// `material` is the index in [`crate::system::System::materials`]. + pub fn new( + label: u32, + index: usize, + kind: ElementKind, + nodes: Vec, + xyz: &[[f64; 3]], + material: usize, + section: RodSection, + ) -> Result { + let RodSection { young, density, area } = section; + let at = |e| Error::element(label, e); + let t = Truss::new(kind, xyz).map_err(at)?; + let omega2 = t.frequency_bound(young, density).map_err(at)?; + let wave_speed = (young / density).sqrt(); + Ok(Rod { + label, + index, + nodes, + material, + young, + stiffness: t.stiffness(young, area).map_err(at)?, + strain: t.strain_operator().map_err(at)?, + nodal_volume: t.lumped_mass(1.0, area).map_err(at)?, + length: 2.0 * wave_speed / omega2.sqrt(), + axial: young * area, + rest_length: t.length().map_err(at)?, + }) + } + + /// The element's degrees of freedom, node by node. + pub fn dofs(&self) -> impl Iterator + '_ { + self.nodes.iter().flat_map(|n| modelcraft_fem::truss::DOFS.node.iter().filter_map(DofKind::direction).map(|d| dof(*n, d))) + } + + /// Add its internal force into `force` and return the element's mean stress and strain + /// tensors: corotational for a rod2 with nonlinear geometry (`xyz` the rest positions of the + /// system's nodes), K u otherwise. + pub fn forces(&self, xyz: &[[f64; 3]], u: &[f64], nonlinear: bool, force: &mut [f64]) -> ([f64; 6], [f64; 6]) { + if nonlinear && let [a, b] = self.nodes[..] { + let at = |n: usize| -> [f64; 3] { + let x = xyz.get(n).copied().unwrap_or_default(); + [0, 1, 2].map(|k| x[k] + u.get(dof(n, k)).copied().unwrap_or(0.0)) + }; + let (pa, pb) = (at(a), at(b)); + let d = [pb[0] - pa[0], pb[1] - pa[1], pb[2] - pa[2]]; + let l = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt(); + let t = if l > 0.0 { d.map(|c| c / l) } else { [0.0; 3] }; + let e = (l - self.rest_length) / self.rest_length; + let n = self.axial * e; + for k in 0..3 { + if let Some(x) = force.get_mut(dof(a, k)) { + *x -= n * t[k]; + } + if let Some(x) = force.get_mut(dof(b, k)) { + *x += n * t[k]; + } + } + return (uniaxial(self.young * e, t, 1.0), uniaxial(e, t, 2.0)); + } + // At most three nodes of three translations: the element vectors live on the stack. + let n = self.stiffness.n.min(MAX_DOFS); + let mut ue = [0.0; MAX_DOFS]; + for (x, d) in ue.iter_mut().zip(self.dofs()) { + *x = u.get(d).copied().unwrap_or(0.0); + } + for (i, d) in self.dofs().enumerate().take(n) { + let f: f64 = (0..n).map(|j| self.stiffness.get(i, j) * ue[j]).sum(); + if let Some(x) = force.get_mut(d) { + *x += f; + } + } + let (mut stress, mut strain) = ([0.0; 6], [0.0; 6]); + let share = 1.0 / self.strain.len().max(1) as f64; + for (row, t) in &self.strain { + let e: f64 = row.iter().zip(&ue).map(|(a, b)| a * b).sum(); + for (k, (s, x)) in uniaxial(self.young * e, *t, 1.0).iter().zip(uniaxial(e, *t, 2.0)).enumerate() { + stress[k] += share * s; + strain[k] += share * x; + } + } + (stress, strain) + } +} + +/// Degrees of freedom of the largest truss (rod3: three nodes of three translations). +const MAX_DOFS: usize = 9; + +/// The global degree of freedom of translation `d` (0–2) of node `n` in the system's uniform, +/// node-major layout ([`NODE_DOFS`] per node, checked when the system is built). +pub fn dof(n: usize, d: usize) -> usize { + NODE_DOFS * n + d +} + +/// The node and the direction (0–2) of a degree of freedom: the inverse of [`dof`]. +pub fn node_direction(dof: usize) -> (usize, usize) { + (dof / NODE_DOFS, dof % NODE_DOFS) +} + +/// A spring or dashpot on degrees of freedom: the elongation is Σ c u(second) − Σ c u(first). +#[derive(Debug, Clone)] +pub struct Link { + pub name: String, + pub dashpot: bool, + /// (degree of freedom, coefficient) with the first end's coefficients negated. + pub terms: Vec<(usize, f64)>, + pub value: f64, + /// For one along the line of its ends: the two nodes and their rest distance. + pub line: Option<(usize, usize, f64)>, +} + +impl Link { + fn stretch(&self, x: &[f64]) -> f64 { + self.terms.iter().map(|(d, c)| c * x.get(*d).copied().unwrap_or(0.0)).sum() + } + + /// Add its internal force (spring on `u`, dashpot on `v`; `xyz` the rest positions of the + /// system's nodes, for a line that turns under nonlinear geometry). + pub fn forces(&self, xyz: &[[f64; 3]], u: &[f64], v: &[f64], nonlinear: bool, stiffness: &mut [f64], viscous: &mut [f64]) { + if nonlinear && let Some((a, b, rest)) = self.line { + let at = |n: usize, x: &[f64]| -> [f64; 3] { [0, 1, 2].map(|k| x.get(dof(n, k)).copied().unwrap_or(0.0)) }; + let pos = |n: usize| { + let (p, d) = (xyz.get(n).copied().unwrap_or_default(), at(n, u)); + [p[0] + d[0], p[1] + d[1], p[2] + d[2]] + }; + let (pa, pb) = (pos(a), pos(b)); + let d = [pb[0] - pa[0], pb[1] - pa[1], pb[2] - pa[2]]; + let l = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt(); + let t = if l > 0.0 { d.map(|c| c / l) } else { [0.0; 3] }; + let (va, vb) = (at(a, v), at(b, v)); + let s = if self.dashpot { self.value * (0..3).map(|k| t[k] * (vb[k] - va[k])).sum::() } else { self.value * (l - rest) }; + let out = if self.dashpot { viscous } else { stiffness }; + for k in 0..3 { + if let Some(x) = out.get_mut(dof(a, k)) { + *x -= s * t[k]; + } + if let Some(x) = out.get_mut(dof(b, k)) { + *x += s * t[k]; + } + } + return; + } + let (x, out) = if self.dashpot { (v, viscous) } else { (u, stiffness) }; + let s = self.value * self.stretch(x); + for (d, c) in &self.terms { + if let Some(f) = out.get_mut(*d) { + *f += s * c; + } + } + } +} + +/// Springs, dashpots and point masses of a model on the explicit system's nodes. +#[derive(Debug, Clone, Default)] +pub struct Discrete { + pub links: Vec, + /// Point mass per degree of freedom. + pub mass: Vec, +} + +impl Discrete { + /// `node` maps a label to its node index. + /// `node` maps a label to its node index, `xyz` gives the nodes' rest positions. + pub fn new(resolved: &Resolved, dofs: usize, xyz: &[[f64; 3]], node: impl Fn(u32) -> Option) -> Result { + let off = |what: &str, name: &str, label: u32| { + Error::NotSupported(format!( + "{what} `{name}` at node {label}: explicit dynamics takes discrete elements on mesh nodes (not reference points) so far" + )) + }; + let mut out = Discrete { links: Vec::new(), mass: vec![0.0; dofs] }; + for m in &resolved.masses { + if m.inertia.iter().any(|c| *c != 0.0) { + return Err(Error::NotSupported(format!("rotary inertia (point mass `{}`) in explicit dynamics: it has no rotations yet", m.name))); + } + let n = node(m.node).ok_or_else(|| off("point mass", &m.name, m.node))?; + for d in 0..3 { + if let Some(x) = out.mass.get_mut(dof(n, d)) { + *x += m.mass[d]; + } + } + } + for s in &resolved.springs { + let what = s.kind(); + if s.direction.is_rotational() { + return Err(Error::NotSupported(format!("rotational {what}s (`{}`) in explicit dynamics: it has no rotations yet", s.name))); + } + let c = s.direction.coefficients(); + let mut terms = Vec::new(); + let mut ends = Vec::new(); + for (label, sign) in std::iter::once((s.first, -1.0)).chain(s.second.map(|b| (b, 1.0))) { + let n = node(label).ok_or_else(|| off(what, &s.name, label))?; + ends.push(n); + terms.extend((0..3).filter(|d| c[*d] != 0.0).map(|d| (dof(n, d), sign * c[d]))); + } + let line = match (s.direction, ends.as_slice()) { + (ResolvedDirection::Line(_), [a, b]) => { + let (pa, pb) = (xyz.get(*a).copied().unwrap_or_default(), xyz.get(*b).copied().unwrap_or_default()); + Some((*a, *b, ((pb[0] - pa[0]).powi(2) + (pb[1] - pa[1]).powi(2) + (pb[2] - pa[2]).powi(2)).sqrt())) + } + _ => None, + }; + out.links.push(Link { name: s.name.clone(), dashpot: s.dashpot, terms, value: s.value, line }); + } + Ok(out) + } + + /// The smallest stable increment of the springs and dashpots on the lumped mass `mass` (per + /// degree of freedom), and the name of a link at the degree of freedom that sets it. A link + /// along its turning line (`nonlinear`) is bounded for every direction it may take. + pub fn critical_increment(&self, mass: &[f64], nonlinear: bool) -> Option<(f64, String)> { + let mut stiff = vec![0.0; mass.len()]; + let mut damp = vec![0.0; mass.len()]; + let mut owner: Vec> = vec![None; mass.len()]; + for (k, l) in self.links.iter().enumerate() { + let rows = if l.dashpot { &mut damp } else { &mut stiff }; + let turning: Vec<(usize, f64)>; + let terms = match l.line { + Some((a, b, _)) if nonlinear => { + turning = (0..3).map(|k| (dof(a, k), -1.0)).chain((0..3).map(|k| (dof(b, k), 1.0))).collect(); + &turning + } + _ => &l.terms, + }; + for (i, ci) in terms { + let mi = mass.get(*i).copied().unwrap_or(0.0); + let sum: f64 = terms.iter().map(|(j, cj)| (l.value * ci * cj).abs() / (mi * mass.get(*j).copied().unwrap_or(0.0)).sqrt()).sum(); + if let Some(r) = rows.get_mut(*i) { + *r += sum; + } + if let Some(o) = owner.get_mut(*i) { + o.get_or_insert(k); + } + } + } + let mut best: Option<(f64, usize)> = None; + for (i, (w2, d)) in stiff.iter().zip(&damp).enumerate() { + let dt = match (*w2 > 0.0, *d > 0.0) { + (false, false) => continue, + (false, true) => 2.0 / d, + (true, _) => { + let w = w2.sqrt(); + let xi = d / (2.0 * w); + (2.0 / w) * ((1.0 + xi * xi).sqrt() - xi) + } + }; + if best.is_none_or(|(b, _)| dt < b) { + best = Some((dt, i)); + } + } + let (dt, i) = best?; + let name = owner.get(i).copied().flatten().and_then(|k| self.links.get(k)).map(|l| l.name.clone()).unwrap_or_default(); + Some((dt, name)) + } +} diff --git a/crates/explicit/src/elements.rs b/crates/explicit/src/elements.rs index 211eb2e..c2e3b2f 100644 --- a/crates/explicit/src/elements.rs +++ b/crates/explicit/src/elements.rs @@ -477,6 +477,22 @@ pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Wor critical = critical.min(element.finish(ws, &es.point.stress, &inc.v, m.bulk)); } } + for r in &sys.rods { + let (sigma, strain) = r.forces(&sys.xyz, m.u, sys.nonlinear, &mut ws.stress_force); + if let Some(es) = state.elements.get_mut(r.index) { + es.point.stress = sigma; + es.strain = strain; + } + if let Some(mat) = sys.materials.get(r.material) { + critical = critical.min(r.length / mat.wave_speed); + } + } + for l in &sys.discrete.links { + l.forces(&sys.xyz, m.u, m.v, sys.nonlinear, &mut ws.stress_force, &mut ws.viscous_force); + } + if let Some((dt, _)) = &sys.discrete_increment { + critical = critical.min(*dt); + } Ok(critical) } diff --git a/crates/explicit/src/integrate.rs b/crates/explicit/src/integrate.rs index 9c78e50..8f9a9b0 100644 --- a/crates/explicit/src/integrate.rs +++ b/crates/explicit/src/integrate.rs @@ -133,13 +133,31 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res .iter() .map(|h| (h.label, h.length, h.material)) .chain(sys.tets.iter().map(|t| (t.label, t.length, t.material))) + .chain(sys.rods.iter().map(|r| (r.label, r.length, r.material))) .map(|(label, length, m)| ControllingElement { label, increment: length / sys.materials.get(m).map_or(f64::INFINITY, |m| m.wave_speed) }) .collect(); critical.sort_by(|a, b| a.increment.total_cmp(&b.increment).then(a.label.cmp(&b.label))); critical.truncate(CONTROLLING); - let smallest = critical - .first() - .map(|c| c.increment) + let n = sys.dofs(); + let mass = sys.dof_mass(); + if let Some(d) = mass.iter().position(|m| !(*m > 0.0 && m.is_finite())) { + let label = sys.labels.get(crate::discrete::node_direction(d).0).copied().unwrap_or(0); + return Err(Error::Invalid(format!( + "node {label} has no mass in direction {}; explicit dynamics needs a mass at every node (give the nodes only springs or dashpots use a point mass)", + crate::discrete::node_direction(d).1 + 1 + ))); + } + // Springs and dashpots have their own stable increment on the lumped mass. + let elements_smallest = critical.first().map(|c| c.increment); + if let Some((l, name)) = &sys.discrete_increment + && elements_smallest.is_none_or(|e| *l < e) + { + warnings.push(format!("spring or dashpot `{name}` limits the time increment to {l:.3e} (before the scale factor)")); + } + let smallest = elements_smallest + .into_iter() + .chain(sys.discrete_increment.as_ref().map(|l| l.0)) + .reduce(f64::min) .filter(|d| *d > 0.0 && d.is_finite()) .ok_or_else(|| Error::Invalid("no element has a finite, positive critical time increment (check the material constants)".into()))?; let estimate = end / (scale * smallest); @@ -170,12 +188,8 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res ))); } - let n = sys.dofs(); - let mass: Vec = (0..n).map(|d| sys.mass.get(d / NODE_DOFS).copied().unwrap_or(0.0)).collect(); - if mass.iter().any(|m| !(*m > 0.0 && m.is_finite())) { - return Err(Error::Internal("a node has no positive mass".into())); - } - let total_mass: f64 = sys.mass.iter().sum(); + // Point masses count with their mean over the three directions. + 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 state = State::new(&sys); @@ -211,17 +225,17 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res history.time.push(t); history.increment.push(dt); history.energies.push(e); - let sum = |x: &[f64], weight: &dyn Fn(usize) -> f64| { + let sum = |x: &[f64], weight: &dyn Fn(usize, usize) -> f64| { let mut s = [0.0; 3]; for (node, c) in x.as_chunks::().0.iter().enumerate() { for k in 0..3 { - s[k] += weight(node) * c[k]; + s[k] += weight(node, k) * c[k]; } } s }; - history.momentum.push(sum(v, &|node| sys.mass.get(node).copied().unwrap_or(0.0))); - history.constraint_force.push(sum(r, &|_| 1.0)); + history.momentum.push(sum(v, &|node, k| mass.get(crate::discrete::dof(node, k)).copied().unwrap_or(0.0))); + history.constraint_force.push(sum(r, &|_, _| 1.0)); let pick = |x: &[f64]| sys.history_nodes.iter().map(|i| [0, 1, 2].map(|k| x.get(NODE_DOFS * i + k).copied().unwrap_or(0.0))).collect(); history.displacement.push(pick(u)); history.velocity.push(pick(v)); diff --git a/crates/explicit/src/lib.rs b/crates/explicit/src/lib.rs index c6dbfb0..75e9f8d 100644 --- a/crates/explicit/src/lib.rs +++ b/crates/explicit/src/lib.rs @@ -32,6 +32,7 @@ // Small fixed-size vector kernels read as the textbook formulas with index loops. #![allow(clippy::needless_range_loop)] +pub mod discrete; pub mod elements; pub mod integrate; pub mod result; diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs index e326e4f..b2b5428 100644 --- a/crates/explicit/src/system.rs +++ b/crates/explicit/src/system.rs @@ -1,5 +1,6 @@ -//! From the resolved model to the explicit system: elements with their materials and stable -//! lengths, lumped masses, prescribed velocities, rigid walls, initial velocities and the load +//! From the resolved model to the explicit system: elements (solids and trusses) with their +//! materials and stable lengths, lumped masses with the point masses, springs and dashpots +//! ([`crate::discrete`]), prescribed velocities, rigid walls, initial velocities and the load //! vector of an explicit dynamic step. use modelcraft_fem::material::Elasticity; @@ -10,6 +11,7 @@ use modelcraft_model::{ Elastic, ElementKind, ExplicitDynamic, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step, StressRate, Target, }; +use crate::discrete::{Discrete, Rod}; use crate::elements::{HexScales, add_at, frequency_bound, hex_scales}; use crate::{Error, Result}; @@ -97,8 +99,14 @@ pub struct System { pub materials: Vec, pub hexes: Vec, pub tets: Vec, + pub rods: Vec, /// Row-sum lumped mass per node. pub mass: Vec, + /// Point masses, springs and dashpots. + pub discrete: Discrete, + /// The springs' and dashpots' stable increment on the lumped mass, and a link that sets it, + /// when there are any. + pub discrete_increment: Option<(f64, String)>, /// Prescribed velocity per degree of freedom (`None` = free). pub prescribed: Vec>, pub walls: Vec, @@ -162,6 +170,7 @@ fn check_formulation(kind: ElementKind, section_name: &str, f: &modelcraft_model ElementKind::Hex8 => Ok(()), ElementKind::Tet4 if f.integration == Integration::Reduced => refuse("reduced integration of tet4 elements"), ElementKind::Tet4 => Ok(()), + ElementKind::Rod2 | ElementKind::Rod3 => Ok(()), other => refuse(&format!("{} elements in explicit dynamics (it has one-point hex8 and nodal-pressure tet4 so far)", other.name())), } } @@ -173,14 +182,18 @@ impl System { } // Elements by label, with their section's material. let mut materials: Vec = Vec::new(); - let mut chosen: Vec<(u32, ElementKind, Vec, usize)> = Vec::new(); + let mut chosen: Vec<(u32, ElementKind, Vec, usize, Option)> = Vec::new(); for (name, section) in model.sections() { - if !matches!(section.kind, SectionKind::Solid {}) { - return Err(Error::NotSupported(format!( - "{} sections (section `{name}`); explicit dynamics has solid elements so far", - section.kind.name() - ))); - } + let area = match section.kind { + SectionKind::Solid {} => None, + SectionKind::Rod { area } => Some(area), + _ => { + return Err(Error::NotSupported(format!( + "{} sections (section `{name}`); explicit dynamics has solid and truss elements so far", + section.kind.name() + ))); + } + }; let mi = match materials.iter().position(|m| m.name == section.material) { Some(i) => i, None => { @@ -191,6 +204,13 @@ impl System { let set = resolved.sets.get(§ion.set).ok_or_else(|| invalid(format!("set `{}` is not resolved", section.set)))?; for &label in &set.elements { let el = model.element(label).ok_or_else(|| invalid(format!("element {label} is missing")))?; + if area.is_some() != matches!(el.kind, ElementKind::Rod2 | ElementKind::Rod3) { + return Err(invalid(format!( + "section `{name}`: a {} section cannot hold {} element {label}", + section.kind.name(), + el.kind.name() + ))); + } check_formulation(el.kind, name, §ion.formulation)?; if el.kind == ElementKind::Tet4 && materials.get(mi).is_some_and(|m| m.j2.is_none()) { return Err(Error::NotSupported(format!( @@ -198,11 +218,13 @@ impl System { section.material ))); } - chosen.push((label, el.kind, el.nodes.clone(), mi)); + chosen.push((label, el.kind, el.nodes.clone(), mi, area)); } } chosen.sort_by_key(|e| e.0); let mut labels: Vec = chosen.iter().flat_map(|e| e.2.iter().copied()).collect(); + // Mesh nodes only discrete elements use (reference points are refused by Discrete::new). + labels.extend(resolved.discrete_labels().filter(|l| model.node(*l).is_some())); labels.sort_unstable(); labels.dedup(); let layout = DofLayout::uniform(labels.iter().copied(), DofSet::TRANSLATIONS); @@ -213,7 +235,7 @@ impl System { return Err(Error::Internal("explicit dynamics needs three translations per node".into())); } if labels.is_empty() { - return Err(invalid("the model has no elements to analyse")); + return Err(invalid("the model has no elements or point masses to analyse")); } let xyz: Vec<[f64; 3]> = labels.iter().map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| invalid(format!("node {l} is missing")))).collect::>()?; @@ -224,7 +246,10 @@ impl System { materials, hexes: Vec::new(), tets: Vec::new(), + rods: Vec::new(), mass: vec![0.0; n], + discrete: Discrete::default(), + discrete_increment: None, prescribed: vec![None; dofs], walls: Vec::new(), initial_velocity: vec![0.0; dofs], @@ -238,7 +263,7 @@ impl System { labels, xyz, }; - for (index, (label, kind, nodes, mi)) in chosen.into_iter().enumerate() { + 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()))?; let (rho, c) = (mat.density, mat.wave_speed); @@ -264,9 +289,34 @@ impl System { sys.lump(&nodes, &mass); sys.tets.push(Tet { label, index, nodes, material: mi, kinematics: k, mass, omega2, length: 2.0 * c / omega2.sqrt() }); } + ElementKind::Rod2 | ElementKind::Rod3 => { + let area = area.ok_or_else(|| Error::Internal("truss area".into()))?; + let source = model.material(&mat.name).ok_or_else(|| Error::Internal("truss material".into()))?; + if source.plastic.is_some() || source.damage.is_some() { + return Err(Error::NotSupported(format!( + "plastic or damaging materials on trusses (material `{}`, element {label}): explicit trusses are elastic so far", + mat.name + ))); + } + if kind == ElementKind::Rod3 && settings.nonlinear_geometry { + return Err(Error::NotSupported(format!( + "three-node trusses with nonlinear geometry (element {label}); set the step's nonlinear_geometry to false" + ))); + } + let young = + modelcraft_fem::truss::young_modulus(source.elastic.as_ref().ok_or_else(|| Error::Internal("truss material".into()))?) + .map_err(|e| invalid(format!("material `{}` of element {label}: {e}", mat.name)))?; + let x: Vec<[f64; 3]> = idx.iter().map(|i| sys.xyz.get(*i).copied().unwrap_or_default()).collect(); + let rod = Rod::new(label, index, kind, idx.clone(), &x, mi, crate::discrete::RodSection { young, density: rho, area })?; + let mass: Vec = rod.nodal_volume.iter().map(|v| rho * v).collect(); + sys.lump(&idx, &mass); + sys.rods.push(rod); + } other => return Err(Error::Internal(format!("{} element kind", other.name()))), } } + sys.discrete = Discrete::new(resolved, sys.dofs(), &sys.xyz, |l| sys.node(l))?; + sys.discrete_increment = sys.discrete.critical_increment(&sys.dof_mass(), sys.nonlinear); sys.prescribe(model, resolved, step, settings.end_time)?; sys.load(model, resolved, step)?; sys.conditions(resolved, settings)?; @@ -290,6 +340,13 @@ impl System { self.layout.dofs() } + /// The lumped mass per degree of freedom: the elements' and the point masses'. + pub fn dof_mass(&self) -> Vec { + (0..self.dofs()) + .map(|d| self.mass.get(crate::discrete::node_direction(d).0).copied().unwrap_or(0.0) + self.discrete.mass.get(d).copied().unwrap_or(0.0)) + .collect() + } + /// Prescribed velocities from the step's boundary conditions: zero for held degrees of /// freedom, value / end time for a displacement reached over the step. fn prescribe(&mut self, model: &Model, resolved: &Resolved, step: &Step, end_time: f64) -> Result<()> { @@ -363,6 +420,12 @@ impl System { add_at(&mut f, *node, &acceleration.map(|a| rho * w * a)); } }); + // Point masses take gravity on the whole model. + if set.is_none() { + for (d, (x, m)) in f.iter_mut().zip(&self.discrete.mass).enumerate() { + *x += m * acceleration.get(crate::discrete::node_direction(d).1).copied().unwrap_or(0.0); + } + } } Load::BodyForce { set, vector } => { let chosen = resolved.sets.get(set).map(|r| r.elements.as_slice()).unwrap_or_default(); @@ -399,6 +462,9 @@ impl System { for t in self.tets.iter().filter(|t| wanted(t.label)) { visit(&t.nodes, &t.kinematics.nodal_volume, density(t.material)); } + for r in self.rods.iter().filter(|r| wanted(r.label)) { + visit(&r.nodes, &r.nodal_volume, density(r.material)); + } } /// The analysed nodes of a resolved set. diff --git a/crates/explicit/tests/one_d.rs b/crates/explicit/tests/one_d.rs new file mode 100644 index 0000000..01ae8b8 --- /dev/null +++ b/crates/explicit/tests/one_d.rs @@ -0,0 +1,307 @@ +//! Trusses, point masses, springs and dashpots in explicit dynamics against closed-form +//! solutions. Run with `--nocapture` to see the errors. +//! +//! - A mass on a spring started with a velocity: x = (v₀/ω) sin ωt. +//! - The same with a dashpot: the amplitude decays as e^{−ξωt}, the period is 2π/ω_d with +//! ω_d = ω√(1 − ξ²); a mass on a dashpot alone slows as v₀ e^{−ct/m}. +//! - Two masses on a two-point spring along its line: ω² = k(1/m₁ + 1/m₂), momentum conserved. +//! - A truss bar driven at one end: a wave front at c t with stress −ρ c v behind it. + +use std::f64::consts::PI; + +use modelcraft_explicit::{ExplicitResult, Limits, solve_explicit}; +use modelcraft_fem::truss::axial; +use modelcraft_model::*; + +fn model(step: ExplicitDynamic) -> Model { + let mut m = Model::default(); + m.add_step(Step::new("Dynamic", Procedure::ExplicitDynamic(ExplicitDynamic::new(step.end_time))), None).unwrap(); + // The step's settings name sets; they are put in once the sets exist (see `solve`). + PENDING.with(|p| *p.borrow_mut() = Some(step)); + m +} + +thread_local! { + static PENDING: std::cell::RefCell> = const { std::cell::RefCell::new(None) }; +} + +/// Put the step's settings in, keeping what is active in it, and solve. +fn solve(m: &mut Model) -> ExplicitResult { + if let Some(settings) = PENDING.with(|p| p.borrow_mut().take()) { + let mut step = m.step("Dynamic").unwrap().clone(); + step.procedure = Procedure::ExplicitDynamic(settings); + m.replace_step(step).unwrap(); + } + solve_explicit(m, None, &Limits::default()).unwrap() +} + +fn set(m: &mut Model, name: &str, nodes: &[u32]) -> Target { + m.add_set(name, Set { nodes: nodes.to_vec(), ..Set::default() }).unwrap(); + Target::Set(name.into()) +} + +fn hold(m: &mut Model, name: &str, target: Target, values: [Option; 6]) { + m.add_boundary_condition(name, BoundaryCondition::Displacement { target, values }).unwrap(); + m.set_boundary_condition_active("Dynamic", name, true).unwrap(); +} + +fn step(end: f64, scale: f64, moving: &str, velocity: [f64; 3]) -> ExplicitDynamic { + let mut s = ExplicitDynamic::new(end); + s.dt_scale = scale; + s.history_nodes = Some(moving.into()); + s.history_interval = Some(end / 20_000.0); + s.initial_velocities = vec![InitialVelocity { set: moving.into(), velocity }]; + s.bulk_viscosity = BulkViscosity { linear: 0.0, quadratic: 0.0 }; + s +} + +/// The history of one direction of the first history node. +fn trace(r: &ExplicitResult, k: usize) -> Vec<(f64, f64)> { + r.history.time.iter().zip(&r.history.displacement).map(|(t, d)| (*t, d[0][k])).collect() +} + +/// Local maxima of a sampled signal (parabolic interpolation through the three samples). +fn peaks(x: &[(f64, f64)]) -> Vec<(f64, f64)> { + x.windows(3) + .filter(|w| w[1].1 > w[0].1 && w[1].1 >= w[2].1) + .map(|w| { + let (h, a, b, c) = (w[1].0 - w[0].0, w[0].1, w[1].1, w[2].1); + let d = 0.5 * (a - c) / (a - 2.0 * b + c); + (w[1].0 + d * h, b - 0.25 * (a - c) * d) + }) + .collect() +} + +/// A mass on a spring to the ground (and optionally a dashpot), started at velocity v₀. +fn oscillator(k: f64, m: f64, damping: Option, periods: f64, scale: f64) -> ExplicitResult { + let omega = (k / m).sqrt(); + let v0 = 0.3; + let mut model = model(step(periods * 2.0 * PI / omega, scale, "n", [v0, 0.0, 0.0])); + model.set_node(1, [0.0, 0.0, 0.0]).unwrap(); + let n = set(&mut model, "n", &[1]); + model.add_spring("k", Spring::Spring { first: n.clone(), second: None, direction: Direction::Dof(1), stiffness: k }).unwrap(); + if let Some(c) = damping { + model.add_spring("c", Spring::Dashpot { first: n.clone(), second: None, direction: Direction::Dof(1), coefficient: c }).unwrap(); + } + model.add_mass("m", PointMass { target: n.clone(), mass: Some(Mass::Isotropic(m)), inertia: None }).unwrap(); + hold(&mut model, "yz", n, [None, Some(0.0), Some(0.0), None, None, None]); + solve(&mut model) +} + +#[test] +fn spring_mass_oscillation() { + let (k, m) = (200.0f64, 0.5); + let r = oscillator(k, m, None, 3.0, 0.02); + let omega = (k / m).sqrt(); + let amplitude = 0.3 / omega; + let worst = trace(&r, 0).iter().map(|(t, x)| (x - amplitude * (omega * t).sin()).abs()).fold(0.0, f64::max); + let p = peaks(&trace(&r, 0)); + let period = (p[2].0 - p[0].0) / 2.0; + eprintln!( + "spring–mass: worst error {:.2e} of the amplitude, period {period:.6} vs {:.6}, Δt {:.3e}, energy error {:.1e}", + worst / amplitude, + 2.0 * PI / omega, + r.initial_increment, + r.energy.relative_error(r.initial_kinetic) + ); + assert!(worst < 2e-3 * amplitude, "second-order phase error"); + assert!((period - 2.0 * PI / omega).abs() < 1e-4 * period); + assert!(r.energy.relative_error(r.initial_kinetic) < 1e-3); + assert!(r.warnings.iter().any(|w| w.contains("limits the time increment")), "{:?}", r.warnings); + // The increment is the scale times the spring's 2/ω. + assert!((r.initial_increment - 0.02 * 2.0 / omega).abs() < 1e-12 * r.initial_increment); +} + +#[test] +fn damped_decay_rate() { + let (k, m, xi) = (200.0f64, 0.5, 0.05); + let omega = (k / m).sqrt(); + let c = 2.0 * xi * m * omega; + let r = oscillator(k, m, Some(c), 6.0, 0.01); + let p = peaks(&trace(&r, 0)); + let wd = omega * (1.0 - xi * xi).sqrt(); + // Decay rate from the first and last peaks, period from their spacing. + let (first, last) = (p[0], p[p.len() - 1]); + let rate = (first.1 / last.1).ln() / (last.0 - first.0); + let period = (last.0 - first.0) / (p.len() - 1) as f64; + let rate_err = (rate - xi * omega).abs() / (xi * omega); + let period_err = (period - 2.0 * PI / wd).abs() / (2.0 * PI / wd); + eprintln!("damped oscillator: decay rate {rate:.6} vs ξω {:.6} ({rate_err:.2e}), period error {period_err:.2e}", xi * omega); + assert!(rate_err < 0.01, "{rate} vs {}", xi * omega); + assert!(period_err < 1e-3); + // The dashpot dissipates what the oscillator loses: the balance holds. + assert!(r.energy.viscous > 0.9 * r.initial_kinetic && r.energy.relative_error(r.initial_kinetic) < 0.01, "{:?}", r.energy); + + // A dashpot alone: v = v₀ e^{−ct/m}. + let mut model = model(step(3.0 * m / c, 0.001, "n", [0.3, 0.0, 0.0])); + model.set_node(1, [0.0; 3]).unwrap(); + let n = set(&mut model, "n", &[1]); + model.add_spring("c", Spring::Dashpot { first: n.clone(), second: None, direction: Direction::Vector([2.0, 0.0, 0.0]), coefficient: c }).unwrap(); + model.add_mass("m", PointMass { target: n.clone(), mass: Some(Mass::Isotropic(m)), inertia: None }).unwrap(); + hold(&mut model, "yz", n, [None, Some(0.0), Some(0.0), None, None, None]); + let r = solve(&mut model); + let (t, v) = (r.history.time.last().unwrap(), r.history.velocity.last().unwrap()[0][0]); + let exact = 0.3 * (-c * t / m).exp(); + eprintln!("dashpot alone: v {v:.6e} vs {exact:.6e} at t = {t:.4}"); + assert!((v - exact).abs() < 0.01 * exact); + // Travel: x(∞) = v₀ m / c. + let x = r.history.displacement.last().unwrap()[0][0]; + let travel = 0.3 * m / c * (1.0 - (-c * t / m).exp()); + assert!((x - travel).abs() < 0.01 * travel, "{x} vs {travel}"); +} + +#[test] +fn two_masses_on_a_line_spring() { + let (k, m1, m2) = (50.0f64, 1.0, 3.0); + let omega = (k * (1.0 / m1 + 1.0 / m2)).sqrt(); + let dir = [0.6, 0.0, 0.8]; + let mut model = model(step(4.0 * 2.0 * PI / omega, 0.02, "a", [0.0; 3])); + model.set_node(1, [0.0; 3]).unwrap(); + model.set_node(2, dir.map(|c| 2.0 * c)).unwrap(); + let a = set(&mut model, "a", &[1]); + let b = set(&mut model, "b", &[2]); + model.add_spring("k", Spring::Spring { first: a.clone(), second: Some(b.clone()), direction: Direction::Line, stiffness: k }).unwrap(); + model.add_mass("m1", PointMass { target: a, mass: Some(Mass::Isotropic(m1)), inertia: None }).unwrap(); + model.add_mass("m2", PointMass { target: b, mass: Some(Mass::Isotropic(m2)), inertia: None }).unwrap(); + // Node 1 starts towards node 2, node 2 at rest: the centre of mass drifts, the spring rings. + PENDING.with(|p| p.borrow_mut().as_mut().unwrap().initial_velocities = vec![InitialVelocity { set: "a".into(), velocity: dir.map(|c| 0.5 * c) }]); + let r = solve(&mut model); + // Momentum is conserved exactly (internal forces are equal and opposite). + let p0 = r.history.momentum[0]; + let drift = r.history.momentum.iter().map(|p| (0..3).map(|k| (p[k] - p0[k]).abs()).fold(0.0, f64::max)).fold(0.0, f64::max); + // The relative stretch along the line oscillates at ω. + let x1 = trace(&r, 0); + let along: Vec<(f64, f64)> = r + .history + .time + .iter() + .zip(&r.history.displacement) + .map(|(t, d)| (*t, d[0][0] * dir[0] + d[0][2] * dir[2] - 0.5 * m1 / (m1 + m2) * t)) + .collect(); + let p = peaks(&along); + let period = (p[3].0 - p[0].0) / 3.0; + eprintln!("two masses: period {period:.6} vs {:.6}, momentum drift {drift:.1e}", 2.0 * PI / omega); + assert!((period - 2.0 * PI / omega).abs() < 1e-4 * period); + assert!(drift < 1e-14); + assert!(x1.iter().all(|(_, x)| x.is_finite())); +} + +/// A rod2 bar along x, held across, its end pushed at velocity v: a wave front at c t and +/// σ = −ρ c v behind it, carried at the bar's own critical increment L/c. +#[test] +fn velocity_wave_in_a_truss_bar() { + let (young, rho, area, length, cells) = (1.0e4f64, 1.0, 0.01, 10.0, 200usize); + let c = (young / rho).sqrt(); + let t = 0.6 * length / c; + let v = 0.01; + for kind in [ElementKind::Rod2, ElementKind::Rod3] { + let mut s = ExplicitDynamic::new(t); + s.output_interval = Some(t / 2.0); + // A rod3 is small-displacement only; the rod2 runs corotational (the default). + s.nonlinear_geometry = kind == ElementKind::Rod2; + let mut m = model(s); + let h = length / cells as f64; + for i in 0..=cells { + m.set_node(1 + i as u32, [i as f64 * h, 0.0, 0.0]).unwrap(); + } + for i in 0..cells { + let mut nodes = vec![1 + i as u32, 2 + i as u32]; + if kind == ElementKind::Rod3 { + m.set_node(1000 + i as u32, [(i as f64 + 0.5) * h, 0.0, 0.0]).unwrap(); + nodes.push(1000 + i as u32); + } + m.set_element(1 + i as u32, kind, nodes).unwrap(); + } + m.add_set("bar", Set { elements: (1..=cells as u32).collect(), ..Set::default() }).unwrap(); + m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young, poisson: 0.0 }), density: Some(rho), ..Material::default() }) + .unwrap(); + m.add_section( + "bar", + Section { set: "bar".into(), material: "Mat".into(), kind: SectionKind::Rod { area }, formulation: Formulation::default() }, + ) + .unwrap(); + let all: Vec = m.nodes().keys().copied().collect(); + let line = set(&mut m, "line", &all); + hold(&mut m, "line", line, [None, Some(0.0), Some(0.0), None, None, None]); + let start = set(&mut m, "start", &[1]); + hold(&mut m, "push", start, [Some(v * t), None, None, None, None, None]); + let r = solve(&mut m); + let f = r.frames.last().unwrap(); + // σ behind the front (x < 0.8 c t) and the front where it falls to half. + let stress: Vec<(f64, f64)> = (0..cells).map(|i| ((i as f64 + 0.5) * h, axial(&f.stress[i], [1.0, 0.0, 0.0], false))).collect(); + let behind: Vec = stress.iter().filter(|(x, _)| *x > 0.1 * c * t && *x < 0.8 * c * t).map(|s| s.1).collect(); + let sigma = behind.iter().sum::() / behind.len() as f64; + let front = + stress.windows(2).find(|w| w[0].1 / (-rho * c * v) >= 0.5 && w[1].1 / (-rho * c * v) < 0.5).map(|w| (w[0].0 + w[1].0) / 2.0).unwrap(); + let expected_dt = if kind == ElementKind::Rod2 { 0.9 * h / c } else { r.initial_increment }; + eprintln!( + "{kind:?} bar wave: σ {sigma:.6} vs −ρcv {:.6}, front {front:.4} vs c t {:.4}, Δt {:.3e}, energy error {:.1e}", + -rho * c * v, + c * t, + r.initial_increment, + r.energy.relative_error(r.initial_kinetic) + ); + assert!((sigma + rho * c * v).abs() < 0.01 * rho * c * v); + assert!((front - c * t).abs() < 3.0 * h); + assert!((r.initial_increment - expected_dt).abs() < 1e-12 * expected_dt); + assert!(r.energy.relative_error(r.initial_kinetic) < 0.01, "{:?}", r.energy); + } +} + +#[test] +fn what_explicit_cannot_take_is_named() { + // A node with a spring and no mass. + let mut m = model(ExplicitDynamic::new(1.0)); + m.set_node(1, [0.0; 3]).unwrap(); + let n = set(&mut m, "n", &[1]); + m.add_spring("k", Spring::Spring { first: n.clone(), second: None, direction: Direction::Dof(1), stiffness: 1.0 }).unwrap(); + let e = solve_explicit(&m, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("node 1 has no mass"), "{e}"); + // Rotary inertia, and a rotational spring. + m.add_reference_point("RP", ReferencePoint { xyz: [0.0; 3], label: 0 }).unwrap(); + m.add_mass("j", PointMass { target: Target::ReferencePoint("RP".into()), mass: None, inertia: Some([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) }).unwrap(); + let e = solve_explicit(&m, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("rotary inertia"), "{e}"); +} + +/// A mass whirling on a stiff rod2 tethered to a fixed point (nonlinear geometry): the rod turns +/// with it and keeps its length, carrying the centripetal force m v²/L. A small-displacement rod +/// would let the mass fly off along the tangent. +#[test] +fn whirling_mass_on_a_corotational_rod() { + let (l, mass, v, axial) = (1.0f64, 1.0, 1.0, 1.0e4); + let t = std::f64::consts::FRAC_PI_2 * l / v; + let mut s = step(t, 0.9, "tip", [0.0, v, 0.0]); + s.bulk_viscosity = BulkViscosity { linear: 0.0, quadratic: 0.0 }; + let mut m = model(s); + m.set_node(1, [0.0; 3]).unwrap(); + m.set_node(2, [l, 0.0, 0.0]).unwrap(); + m.set_element(1, ElementKind::Rod2, vec![1, 2]).unwrap(); + m.add_set("rod", Set { elements: vec![1], ..Set::default() }).unwrap(); + // A light rod (E A = 1e4, ρ A = 1e-4) with the mass at its end. + m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: axial, poisson: 0.0 }), density: Some(1e-4), ..Material::default() }) + .unwrap(); + m.add_section( + "rod", + Section { set: "rod".into(), material: "Mat".into(), kind: SectionKind::Rod { area: 1.0 }, formulation: Formulation::default() }, + ) + .unwrap(); + let pivot = set(&mut m, "pivot", &[1]); + m.add_boundary_condition("pivot", BoundaryCondition::Fixed { target: pivot }).unwrap(); + m.set_boundary_condition_active("Dynamic", "pivot", true).unwrap(); + let tip = set(&mut m, "tip", &[2]); + m.add_mass("m", PointMass { target: tip.clone(), mass: Some(Mass::Isotropic(mass)), inertia: None }).unwrap(); + hold(&mut m, "plane", tip, [None, None, Some(0.0), None, None, None]); + let r = solve(&mut m); + let radius = r.history.displacement.iter().map(|d| ((l + d[0][0]).powi(2) + d[0][1].powi(2)).sqrt()).fold(0.0f64, |w, x| w.max((x - l).abs())); + let end = r.history.displacement.last().unwrap()[0]; + eprintln!( + "whirling mass: largest change of radius {radius:.3e} (centripetal stretch m v²/EA = {:.1e}), end at ({:.4}, {:.4})", + mass * v * v / axial, + l + end[0], + end[1] + ); + assert!(radius < 5.0 * mass * v * v / axial, "the rod keeps its length to its elastic stretch"); + // A quarter turn: the mass ends near (0, L). + assert!((l + end[0]).abs() < 1e-2 && (end[1] - l).abs() < 1e-2, "{end:?}"); +} diff --git a/crates/fem/src/basis.rs b/crates/fem/src/basis.rs index 5ee7c47..ff1d6eb 100644 --- a/crates/fem/src/basis.rs +++ b/crates/fem/src/basis.rs @@ -199,6 +199,7 @@ impl Basis { pub fn of(kind: ElementKind) -> Basis { match kind { ElementKind::Rod2 | ElementKind::Beam2 => Basis::Line2, + ElementKind::Rod3 => Basis::Line3, ElementKind::Tri3 => Basis::Tri3, ElementKind::Tri6 => Basis::Tri6, ElementKind::Quad4 => Basis::Quad4, diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs index f9decff..42f8a1d 100644 --- a/crates/fem/src/lib.rs +++ b/crates/fem/src/lib.rs @@ -24,6 +24,8 @@ //! Green–Naghdi); //! - [`plasticity`]: J2 radial return with isotropic, kinematic and Johnson–Cook hardening, //! Johnson–Cook damage; +//! - [`truss`]: two- and three-node trusses (straight or curved), their stiffness, mass, body +//! force and axial stress; //! - [`dense`]: the small dense linear algebra inside an element. //! //! Conventions: element matrices are ordered node by node, each node's degrees of freedom in @@ -46,6 +48,7 @@ pub mod plasticity; pub mod quadrature; pub mod rotation; pub mod shape; +pub mod truss; pub mod uniform; pub use basis::Basis; diff --git a/crates/fem/src/truss.rs b/crates/fem/src/truss.rs new file mode 100644 index 0000000..0c11b60 --- /dev/null +++ b/crates/fem/src/truss.rs @@ -0,0 +1,404 @@ +//! Truss (axial) elements: the two-node rod2 and the three-node rod3, straight or curved, in any +//! orientation (a plane truss is a space truss whose nodes lie in a plane). +//! +//! **Formulation.** The isoparametric bar of R. D. Cook, D. S. Malkus, M. E. Plesha, R. J. Witt, +//! *Concepts and Applications of Finite Element Analysis*, 4th ed., Wiley 2002, §3.6 and §6.2 +//! (and T. J. R. Hughes, *The Finite Element Method*, Dover 2000, §1.10), written along the +//! member's arc in three dimensions. On ξ ∈ [−1, 1] the shape functions are +//! N = ((1 − ξ)/2, (1 + ξ)/2) for rod2 and N = (ξ(ξ − 1)/2, ξ(ξ + 1)/2, 1 − ξ²) for rod3 (the +//! ends, then the middle node). The position is x(ξ) = Σ Nᵢ xᵢ, the arc length element +//! ds = J dξ with J = |x′(ξ)|, and the unit tangent t = x′/J. The axial strain is the stretch +//! along the tangent, ε = t · ∂u/∂s = Σᵢ (Nᵢ′/J) t · uᵢ, so with Bᵢ = (Nᵢ′/J) tᵀ: +//! +//! - stiffness K = ∫ E A Bᵀ B ds, +//! - consistent mass ρ A ∫ Nᵢ Nⱼ ds per direction, lumped mass by HRZ diagonal scaling (E. Hinton, +//! T. Rock, O. C. Zienkiewicz, *A note on mass lumping and related processes in the finite +//! element method*, Earthq. Eng. Struct. Dyn. 4(3):245–249, 1976), which gives ρAL/2 at each +//! end of rod2 and ρAL·(1/6, 1/6, 2/3) on a straight rod3, +//! - body force A ∫ Nᵢ b ds for a force b per volume, +//! - stress σ = E ε along the tangent: in global axes the uniaxial tensor σ t tᵀ (strain ε t tᵀ, +//! with engineering shear strains), so stresses of trusses and solids read the same way and +//! the axial force is N = A tᵀ σ t. +//! +//! The rules are Gauss–Legendre with two points for rod2 and three for rod3: exact for the +//! stiffness and mass of a straight member, and the usual full integration of a curved one. +//! Material: the Young's modulus of an isotropic material; the truss knows no other constant. + +use modelcraft_model::{Elastic, ElementKind}; + +use crate::basis::Basis; +use crate::dense::Matrix; +use crate::dof::{DofSet, ElementDofs}; +use crate::element::PointStress; +use crate::quadrature::gauss_legendre; +use crate::{Error, Result}; + +/// What a truss's nodes carry: the three translations (a plane truss is a space truss whose +/// out-of-plane motion is held; trusses of planar models will declare [`DofSet::PLANE`]). +pub const DOFS: ElementDofs = ElementDofs::new(DofSet::TRANSLATIONS); + +/// The Young's modulus a truss uses, or why there is none. +pub fn young_modulus(elastic: &Elastic) -> Result { + match *elastic { + Elastic::Isotropic { young, .. } if young > 0.0 && young.is_finite() => Ok(young), + Elastic::Isotropic { young, .. } => Err(Error::Invalid(format!("Young's modulus must be positive, got {young}"))), + Elastic::EngineeringConstants { .. } => { + Err(Error::NotSupported("orthotropic materials on trusses (a truss takes one Young's modulus: use an isotropic material)".into())) + } + } +} + +/// Shape functions and their ξ-derivatives ([`Basis::Line2`], [`Basis::Line3`]). +fn shape(n: usize, xi: f64) -> ([f64; 3], [f64; 3]) { + let s = if n == 2 { Basis::Line2 } else { Basis::Line3 }.eval([xi, 0.0, 0.0]); + let (mut v, mut d) = ([0.0; 3], [0.0; 3]); + for i in 0..n.min(3) { + v[i] = s.n[i]; + d[i] = s.dn[i][0]; + } + (v, d) +} + +/// The geometry at one point: shape functions, their arc-length derivatives Nᵢ′/J, the unit +/// tangent and the arc length per unit ξ. +#[derive(Debug, Clone)] +struct Point { + n: [f64; 3], + ds: [f64; 3], + tangent: [f64; 3], + jacobian: f64, + weight: f64, + xyz: [f64; 3], +} + +/// A truss element on its node coordinates (model node order), with the geometry at its +/// integration points worked out once. +#[derive(Debug, Clone)] +pub struct Truss<'a> { + pub kind: ElementKind, + xyz: &'a [[f64; 3]], + rule: Vec, +} + +impl<'a> Truss<'a> { + pub fn new(kind: ElementKind, xyz: &'a [[f64; 3]]) -> Result> { + if !matches!(kind, ElementKind::Rod2 | ElementKind::Rod3) { + return Err(Error::Invalid(format!("{} is not a truss element", kind.name()))); + } + if xyz.len() != kind.node_count() { + return Err(Error::Invalid(format!("{} needs {} nodes, got {}", kind.name(), kind.node_count(), xyz.len()))); + } + // The tangent must not vanish or turn back along the member (a middle node off its span). + let scale = chord(xyz).max(f64::MIN_POSITIVE); + let reference = at(xyz, 0.0, 0.0)?.tangent; + for xi in [-1.0, -0.5, 0.5, 1.0] { + let p = at(xyz, xi, 0.0)?; + if !(p.jacobian > 1e-9 * scale) || p.tangent.iter().zip(&reference).map(|(a, b)| a * b).sum::() <= 0.0 { + return Err(Error::Degenerate("the member folds back on itself or has no length".into())); + } + } + let rule = gauss_legendre(kind.node_count()).iter().map(|(xi, w)| at(xyz, *xi, *w)).collect::>()?; + Ok(Truss { kind, xyz, rule }) + } + + fn nodes(&self) -> usize { + self.kind.node_count() + } + + /// The integration points: two for rod2, three for rod3. + fn points(&self) -> Result<&[Point]> { + Ok(&self.rule) + } + + /// The natural coordinates of the integration points, in order. + pub fn abscissae(&self) -> Vec { + gauss_legendre(self.nodes()).iter().map(|(xi, _)| *xi).collect() + } + + /// The arc length. + pub fn length(&self) -> Result { + Ok(self.points()?.iter().map(|p| p.weight * p.jacobian).sum()) + } + + /// The stiffness (3n × 3n, translations node by node) for Young's modulus `young` and + /// cross-sectional area `area`. + pub fn stiffness(&self, young: f64, area: f64) -> Result { + let n = self.nodes(); + let mut k = Matrix::zeros(3 * n); + for p in self.points()? { + let b: Vec = (0..3 * n).map(|a| p.ds[a / 3] * p.tangent[a % 3]).collect(); + let factor = young * area * p.weight * p.jacobian; + for (i, bi) in b.iter().enumerate() { + for (j, bj) in b.iter().enumerate() { + k.add(i, j, factor * bi * bj); + } + } + } + Ok(k) + } + + /// The consistent mass per direction (n × n). + pub fn mass(&self, density: f64, area: f64) -> Result { + let n = self.nodes(); + let mut m = Matrix::zeros(n); + for p in self.points()? { + let factor = density * area * p.weight * p.jacobian; + for i in 0..n { + for j in 0..n { + m.add(i, j, factor * p.n[i] * p.n[j]); + } + } + } + Ok(m) + } + + /// The HRZ lumped mass per node: the consistent diagonal scaled to the member's mass. + pub fn lumped_mass(&self, density: f64, area: f64) -> Result> { + let m = self.mass(density, area)?; + let n = self.nodes(); + let diagonal: Vec = (0..n).map(|i| m.get(i, i)).collect(); + let total = density * area * self.length()?; + let sum: f64 = diagonal.iter().sum(); + if !(sum > 0.0) { + return Err(Error::Degenerate("the member has no mass to lump".into())); + } + Ok(diagonal.iter().map(|d| d * total / sum).collect()) + } + + /// Consistent nodal forces of a force `b` per volume. + pub fn body_force(&self, b: [f64; 3], area: f64) -> Result> { + let mut f = vec![[0.0; 3]; self.nodes()]; + for p in self.points()? { + for (i, fi) in f.iter_mut().enumerate() { + let w = area * p.weight * p.jacobian * p.n[i]; + for k in 0..3 { + fi[k] += w * b[k]; + } + } + } + Ok(f) + } + + /// At each integration point, the row that gives the axial strain from the element's + /// displacements (3n, node by node), and the tangent there: ε = row · u. + pub fn strain_operator(&self) -> Result, [f64; 3])>> { + let n = self.nodes(); + Ok(self.points()?.iter().map(|p| ((0..3 * n).map(|a| p.ds[a / 3] * p.tangent[a % 3]).collect(), p.tangent)).collect()) + } + + /// The axial strain at each integration point for nodal displacements `u`. + pub fn axial_strains(&self, u: &[[f64; 3]]) -> Result> { + if u.len() != self.nodes() { + return Err(Error::Invalid(format!("{} displacements for a {}", u.len(), self.kind.name()))); + } + Ok(self.points()?.iter().map(|p| (0..self.nodes()).map(|i| p.ds[i] * (0..3).map(|k| p.tangent[k] * u[i][k]).sum::()).sum()).collect()) + } + + /// Strain and stress at the integration points as global tensors (see the module docs). + pub fn stresses(&self, young: f64, u: &[[f64; 3]]) -> Result> { + let strains = self.axial_strains(u)?; + Ok(self + .points()? + .iter() + .zip(strains) + .map(|(p, e)| PointStress { xyz: p.xyz, strain: uniaxial(e, p.tangent, 2.0), stress: uniaxial(young * e, p.tangent, 1.0) }) + .collect()) + } + + /// Integration-point values extrapolated to the nodes (the polynomial through the points, + /// evaluated at the nodes' ξ: −1, 1, then 0). + pub fn extrapolate(&self, values: &[[f64; 6]]) -> Result> { + let xs = self.abscissae(); + if values.len() != xs.len() { + return Err(Error::Invalid(format!("{} values for {} integration points", values.len(), xs.len()))); + } + let at = |xi: f64| { + let mut out = [0.0; 6]; + for (a, va) in xs.iter().zip(values) { + let l: f64 = xs.iter().filter(|b| *b != a).map(|b| (xi - b) / (a - b)).product(); + for k in 0..6 { + out[k] += l * va[k]; + } + } + out + }; + Ok([-1.0, 1.0, 0.0].into_iter().take(self.nodes()).map(at).collect()) + } + + /// The tangent at each integration point. + pub fn tangents(&self) -> Result> { + Ok(self.points()?.iter().map(|p| p.tangent).collect()) + } + + /// An upper bound of the highest natural frequency squared of the element with its lumped + /// mass, for explicit time integration: exact (4c²/L²) for rod2, and Gershgorin's bound on + /// M^{-1/2} K M^{-1/2} for rod3. + pub fn frequency_bound(&self, young: f64, density: f64) -> Result { + if self.kind == ElementKind::Rod2 { + let l = chord(self.xyz); + return Ok(4.0 * young / (density * l * l)); + } + // The area cancels: K ∝ A and M ∝ A. + let k = self.stiffness(young, 1.0)?; + let m = self.lumped_mass(density, 1.0)?; + let mut worst: f64 = 0.0; + for i in 0..k.n { + let mi = m.get(i / 3).copied().unwrap_or(0.0); + let row: f64 = (0..k.n).map(|j| k.get(i, j).abs() / (mi * m.get(j / 3).copied().unwrap_or(0.0)).sqrt()).sum(); + worst = worst.max(row); + } + Ok(worst) + } +} + +/// The distance between the ends. +fn chord(xyz: &[[f64; 3]]) -> f64 { + match (xyz.first(), xyz.get(1)) { + (Some(a), Some(b)) => ((a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2) + (a[2] - b[2]).powi(2)).sqrt(), + _ => 0.0, + } +} + +/// The geometry at `xi` of the element on `xyz` (two or three nodes). +fn at(xyz: &[[f64; 3]], xi: f64, weight: f64) -> Result { + let (n, dn) = shape(xyz.len(), xi); + let mut d = [0.0; 3]; + let mut x = [0.0; 3]; + for (i, p) in xyz.iter().enumerate() { + for k in 0..3 { + d[k] += dn[i] * p[k]; + x[k] += n[i] * p[k]; + } + } + let jacobian = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt(); + if !(jacobian > 0.0 && jacobian.is_finite()) { + return Err(Error::Degenerate("the member has no length".into())); + } + Ok(Point { n, ds: dn.map(|v| v / jacobian), tangent: d.map(|c| c / jacobian), jacobian, weight, xyz: x }) +} + +/// The tensor s t tᵀ as a 6-vector (xx, yy, zz, xy, xz, yz), shear terms times `shear` (2 for +/// engineering shear strains). +pub fn uniaxial(s: f64, t: [f64; 3], shear: f64) -> [f64; 6] { + [s * t[0] * t[0], s * t[1] * t[1], s * t[2] * t[2], shear * s * t[0] * t[1], shear * s * t[0] * t[2], shear * s * t[1] * t[2]] +} + +/// The axial value (stress or strain) of a uniaxial tensor along a unit tangent: tᵀ σ t, with +/// engineering shear strains halved when `engineering`. +pub fn axial(v: &[f64; 6], t: [f64; 3], engineering: bool) -> f64 { + let h = if engineering { 0.5 } else { 1.0 }; + v[0] * t[0] * t[0] + v[1] * t[1] * t[1] + v[2] * t[2] * t[2] + 2.0 * h * (v[3] * t[0] * t[1] + v[4] * t[0] * t[2] + v[5] * t[1] * t[2]) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn norm(a: [f64; 3]) -> f64 { + (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt() + } + + /// A rod along a skew direction, and a rod3 with its middle node at the midpoint. + fn straight(kind: ElementKind) -> Vec<[f64; 3]> { + let (a, b) = ([0.5, -1.0, 2.0], [2.5, 1.0, 3.0]); + let mut x = vec![a, b]; + if kind == ElementKind::Rod3 { + x.push([0, 1, 2].map(|k| (a[k] + b[k]) / 2.0)); + } + x + } + + #[test] + fn straight_rods_match_the_textbook_matrices() { + let (young, area, density) = (200.0, 3.0, 2.0); + for kind in [ElementKind::Rod2, ElementKind::Rod3] { + let x = straight(kind); + let t = Truss::new(kind, &x).unwrap(); + let l = norm([2.0, 2.0, 1.0]); + assert!((t.length().unwrap() - l).abs() < 1e-14); + // Rigid translations and the rotation about the axis' normals: zero force. + let k = t.stiffness(young, area).unwrap(); + for d in 0..3 { + let u: Vec = (0..k.n).map(|a| if a % 3 == d { 1.0 } else { 0.0 }).collect(); + assert!(k.mul_vec(&u).iter().all(|f| f.abs() < 1e-12), "{kind:?}"); + } + // A uniform stretch ε along the axis: strain ε at every point, end forces ±EAε. + let eps = 1e-3; + let dir = [2.0 / l, 2.0 / l, 1.0 / l]; + let u: Vec<[f64; 3]> = + x.iter().map(|p| [0, 1, 2].map(|c| eps * dir[c] * ((p[0] - 0.5) * dir[0] + (p[1] + 1.0) * dir[1] + (p[2] - 2.0) * dir[2]))).collect(); + for e in t.axial_strains(&u).unwrap() { + assert!((e - eps).abs() < 1e-15, "{kind:?} {e}"); + } + let flat: Vec = u.iter().flatten().copied().collect(); + let f = k.mul_vec(&flat); + for c in 0..3 { + assert!((f[3 + c] - young * area * eps * dir[c]).abs() < 1e-12, "{kind:?}"); + } + let s = t.stresses(young, &u).unwrap(); + for p in &s { + assert!((axial(&p.stress, dir, false) - young * eps).abs() < 1e-12); + assert!((axial(&p.strain, dir, true) - eps).abs() < 1e-15); + } + let nodal = t.extrapolate(&s.iter().map(|p| p.stress).collect::>()).unwrap(); + assert!(nodal.iter().all(|v| (axial(v, dir, false) - young * eps).abs() < 1e-12)); + // Mass: total ρAL per direction; lumped HRZ. + let m = t.mass(density, area).unwrap(); + let total = density * area * l; + assert!((m.data.iter().sum::() - total).abs() < 1e-12); + let lumped = t.lumped_mass(density, area).unwrap(); + let want: Vec = if kind == ElementKind::Rod2 { vec![0.5, 0.5] } else { vec![1.0 / 6.0, 1.0 / 6.0, 2.0 / 3.0] }; + for (a, b) in lumped.iter().zip(want) { + assert!((a - b * total).abs() < 1e-12, "{kind:?} {lumped:?}"); + } + let fb = t.body_force([0.0, 0.0, -1.0], area).unwrap(); + assert!((fb.iter().map(|f| f[2]).sum::() + area * l).abs() < 1e-12); + // The highest frequency of the free element with its lumped mass: exact for rod2, + // an upper bound of the exact one for rod3. + let bound = t.frequency_bound(young, density).unwrap(); + let c2 = young / density; + if kind == ElementKind::Rod2 { + assert!((bound - 4.0 * c2 / (l * l)).abs() < 1e-9 * bound); + } else { + // Exact: axial eigenproblem of K = EA/(3L)[7 1 −8; 1 7 −8; −8 −8 16] with + // M = ρAL diag(1/6, 1/6, 2/3): ω²max = 24 c²/L² (antisymmetric 12, symmetric 24). + assert!(bound >= 24.0 * c2 / (l * l) * (1.0 - 1e-12), "{bound} vs {}", 24.0 * c2 / (l * l)); + } + } + } + + #[test] + fn bad_geometry_is_an_error() { + assert!(Truss::new(ElementKind::Rod2, &[[1.0; 3], [1.0; 3]]).is_err()); + // The middle node beyond an end folds the member back. + assert!(Truss::new(ElementKind::Rod3, &[[0.0; 3], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]]).is_err()); + assert!(Truss::new(ElementKind::Rod3, &[[0.0; 3], [1.0, 0.0, 0.0]]).is_err()); + assert!(Truss::new(ElementKind::Hex8, &[[0.0; 3]; 8]).is_err()); + assert!(Truss::new(ElementKind::Rod2, &[[0.0; 3], [f64::NAN, 0.0, 0.0]]).is_err()); + assert!(young_modulus(&Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [1.0; 3] }).is_err()); + } + + /// A curved member: a circular arc of radius R (middle node on the arc). Any linear + /// displacement field is reproduced exactly: a uniform expansion u = s x (a ring growing + /// radially by δ = s R) strains it by s everywhere, an infinitesimal rotation not at all. + #[test] + fn curved_rod3_reproduces_linear_fields() { + let (radius, young) = (3.0, 70.0); + for span in [0.3f64, 1.0, 2.0] { + let at = |a: f64| [radius * a.cos(), radius * a.sin(), 0.5]; + let x = [at(0.0), at(span), at(span / 2.0)]; + let t = Truss::new(ElementKind::Rod3, &x).unwrap(); + let (s, w) = (2e-3, [1e-3, -2e-3, 3e-3]); + let u: Vec<[f64; 3]> = + x.iter().map(|p| [0, 1, 2].map(|k| s * p[k] + w[(k + 1) % 3] * p[(k + 2) % 3] - w[(k + 2) % 3] * p[(k + 1) % 3])).collect(); + let points = t.stresses(young, &u).unwrap(); + for (p, tangent) in points.iter().zip(t.tangents().unwrap()) { + assert!((axial(&p.strain, tangent, true) - s).abs() < 1e-15, "{span}"); + assert!((axial(&p.stress, tangent, false) - young * s).abs() < 1e-13, "{span}"); + } + // The arc length: exact for a quadratic arc only as the span shrinks. + let err = (t.length().unwrap() - radius * span).abs() / (radius * span); + assert!(err < 0.01 * span.powi(4), "{span}: {err}"); + } + } +} diff --git a/crates/io/src/mesh.rs b/crates/io/src/mesh.rs index f21a285..8a78849 100644 --- a/crates/io/src/mesh.rs +++ b/crates/io/src/mesh.rs @@ -38,7 +38,7 @@ const PART: &str = "mesh.bin"; const NO_TAG: u32 = u32::MAX; /// The element kind codes of the format. -pub const KIND_CODES: [(ElementKind, u8); 11] = [ +pub const KIND_CODES: [(ElementKind, u8); 12] = [ (ElementKind::Rod2, 1), (ElementKind::Beam2, 2), (ElementKind::Tri3, 3), @@ -50,6 +50,7 @@ pub const KIND_CODES: [(ElementKind, u8); 11] = [ (ElementKind::Wedge6, 9), (ElementKind::Hex8, 10), (ElementKind::Hex20, 11), + (ElementKind::Rod3, 12), ]; fn kind_code(kind: ElementKind) -> u8 { diff --git a/crates/mesh/src/quality/mod.rs b/crates/mesh/src/quality/mod.rs index 18ceecf..efc99e3 100644 --- a/crates/mesh/src/quality/mod.rs +++ b/crates/mesh/src/quality/mod.rs @@ -205,7 +205,9 @@ pub fn element_quality(kind: ElementKind, xyz: &[[f64; 3]]) -> Result return Err(Error::Invalid(format!("quality metrics are not defined for {}", kind.name()))), + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => { + return Err(Error::Invalid(format!("quality metrics are not defined for {}", kind.name()))); + } ElementKind::Tet4 | ElementKind::Tet10 => shape::tet(corners(&p)?), ElementKind::Wedge6 => { let c: [[f64; 3]; 6] = corners(&p)?; diff --git a/crates/model/src/check.rs b/crates/model/src/check.rs index 4146575..76b881b 100644 --- a/crates/model/src/check.rs +++ b/crates/model/src/check.rs @@ -71,6 +71,11 @@ pub struct Resolved { /// apply, so every one of them distributes a load the same way. A force or traction gives /// force components, a moment moment components. pub nodal_loads: BTreeMap>, + /// Every point mass at every node of its target, by definition name, then node. + pub masses: Vec, + /// Every spring and dashpot at every node of its target (or between its two ends), by + /// definition name, then node. + pub springs: Vec, /// The check's warnings (a resolved model has no errors). pub warnings: Vec, /// How each element of a 2D section (plane stress, plane strain, axisymmetric) extends out @@ -97,6 +102,11 @@ impl Resolved { self.sets.get(set).map(|s| s.edges.iter().filter(|e| self.planar.contains_key(&e.element)).copied().collect()).unwrap_or_default() } + /// The node labels point masses, springs and dashpots act on (repeats possible). + pub fn discrete_labels(&self) -> impl Iterator + '_ { + self.masses.iter().map(|m| m.node).chain(self.springs.iter().flat_map(|s| std::iter::once(s.first).chain(s.second))) + } + /// The nodes a boundary condition, force or moment acts on. pub fn target_nodes(&self, target: &Target) -> &[u32] { match target { @@ -160,9 +170,10 @@ impl Model { let mut resolved = Resolved::default(); let units = self.units; - // A model of fluid steps alone needs no mesh: its domains carry their own grid. + // A model of fluid steps alone needs no mesh: its domains carry their own grid; one of + // point masses, springs and dashpots alone is its own mesh. let fluid_only = !self.steps.is_empty() && self.steps.iter().all(|s| s.procedure.fluid_domain().is_some()); - if self.elements.is_empty() && !fluid_only { + if self.elements.is_empty() && !fluid_only && self.masses.is_empty() && self.springs.is_empty() { d.error("model", "the mesh has no elements"); } if units.length().is_none() { @@ -229,7 +240,7 @@ impl Model { k.dimension() == 2 } SectionKind::Beam { .. } => k == ElementKind::Beam2, - SectionKind::Rod { .. } => k == ElementKind::Rod2, + SectionKind::Rod { .. } => matches!(k, ElementKind::Rod2 | ElementKind::Rod3), }; // The linear element with an incompatible-modes formulation for this kind of section. let incompatible_kind = if sec.kind.is_planar() { ElementKind::Quad4 } else { ElementKind::Hex8 }; @@ -362,6 +373,9 @@ impl Model { } } + // Point masses, springs and dashpots. + let (masses, springs) = self.resolve_discrete(&resolved, &mut d, &has_rotations); + // Loads. let mut nodal_loads = BTreeMap::new(); for (name, load) in &self.loads { @@ -449,11 +463,123 @@ impl Model { } resolved.nodal_loads = nodal_loads; + resolved.masses = masses; + resolved.springs = springs; let mut out = d.0; out.sort_by_key(|a| a.severity); (resolved, out) } + /// Point masses, springs and dashpots expanded to their nodes, reporting what cannot be + /// expanded. + fn resolve_discrete( + &self, + resolved: &Resolved, + d: &mut Diagnostics, + has_rotations: &dyn Fn(&Target) -> bool, + ) -> (Vec, Vec) { + use crate::{Direction, ResolvedDirection, ResolvedMass, ResolvedSpring}; + let mut masses = Vec::new(); + let mut springs = Vec::new(); + let mut count = 0usize; + let mut over = false; + let mut take = |n: usize| { + count = count.saturating_add(n); + over |= count > crate::discrete::MAX_DISCRETE; + !over + }; + for (name, m) in &self.masses { + let subject = format!("point mass `{name}`"); + let nodes = resolved.target_nodes(&m.target); + if nodes.is_empty() { + d.error(&subject, format!("{} has no nodes", m.target)); + continue; + } + if m.inertia.is_some_and(|i| i.iter().any(|c| *c != 0.0)) && !has_rotations(&m.target) { + d.error(&subject, format!("{} has no nodes with rotations for its inertia; put it on a reference point", m.target)); + } + if !take(nodes.len()) { + break; + } + let mass = m.mass.map_or([0.0; 3], crate::Mass::components); + let inertia = m.inertia.unwrap_or_default(); + masses.extend(nodes.iter().map(|n| ResolvedMass { name: name.clone(), node: *n, mass, inertia })); + } + for (name, s) in &self.springs { + let subject = format!("{} `{name}`", s.kind()); + let first = resolved.target_nodes(s.first()); + let rotational = s.direction().is_rotational(); + if first.is_empty() { + d.error(&subject, format!("{} has no nodes", s.first())); + continue; + } + if rotational && s.targets().any(|t| !has_rotations(t)) { + d.error(&subject, "acts on a rotation, and its nodes have none; connect reference points (or shell and beam nodes)"); + } + let unit = |v: [f64; 3]| { + let len = (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt(); + (len > 0.0 && len.is_finite()).then(|| v.map(|c| c / len)) + }; + let pairs: Vec<(u32, Option)> = match s.second() { + None => first.iter().map(|n| (*n, None)).collect(), + Some(t) => { + let second = resolved.target_nodes(t); + if first.len() != 1 || second.len() != 1 { + d.error( + &subject, + format!( + "a two-point {} connects one node to one node; {} has {} and {t} has {} (use a vertex set or a reference point at each end)", + s.kind(), + s.first(), + first.len(), + second.len() + ), + ); + continue; + } + match (first.first(), second.first()) { + (Some(a), Some(b)) if a == b => { + d.error(&subject, format!("both ends are node {a}")); + continue; + } + (Some(a), Some(b)) => vec![(*a, Some(*b))], + _ => continue, + } + } + }; + if !take(pairs.len()) { + break; + } + for (a, b) in pairs { + let direction = match s.direction() { + Direction::Dof(k) => ResolvedDirection::Dof(k), + Direction::Vector(v) => match unit(v) { + Some(n) => ResolvedDirection::Vector(n), + None => { + d.error(&subject, "its direction vector has no length"); + break; + } + }, + Direction::Line => { + let ends = b.and_then(|b| Some((self.label_xyz(a)?, self.label_xyz(b)?))); + match ends.and_then(|(x, y)| unit([y[0] - x[0], y[1] - x[1], y[2] - x[2]])) { + Some(n) => ResolvedDirection::Line(n), + None => { + d.error(&subject, "its two ends are at the same point, so it has no line to act along"); + break; + } + } + } + }; + springs.push(ResolvedSpring { name: name.clone(), dashpot: s.is_dashpot(), first: a, second: b, direction, value: s.value() }); + } + } + if over { + d.error("model", format!("point masses, springs and dashpots resolve to more than {} copies", crate::discrete::MAX_DISCRETE)); + } + (masses, springs) + } + /// A force, moment or traction as nodal values (`None` for other loads), or why it cannot be /// distributed. fn nodal_values(&self, load: &Load, r: &Resolved) -> std::result::Result>, String> { @@ -505,13 +631,17 @@ impl Model { })) } - /// The rules of a 2D model: one kind of 2D model (plane or axisymmetric) and nothing else, in - /// the XY plane (x ≥ 0 when axisymmetric), elements numbered counter-clockwise. + /// The rules of a 2D model: one kind of 2D model (plane or axisymmetric) and nothing else but + /// trusses (plane models only), in the XY plane (x ≥ 0 when axisymmetric), elements numbered + /// counter-clockwise. fn check_planar(&self, r: &Resolved, d: &mut Diagnostics) { if !r.is_planar() { return; } - let others: Vec<&str> = self.sections.iter().filter(|(_, s)| !s.kind.is_planar()).map(|(n, _)| n.as_str()).collect(); + let axisymmetric = r.extent() == Some(Extent::Revolved); + // Trusses join a plane model with its two in-plane unknowns per node. + let truss = |s: &crate::Section| matches!(s.kind, SectionKind::Rod { .. }) && !axisymmetric; + let others: Vec<&str> = self.sections.iter().filter(|(_, s)| !s.kind.is_planar() && !truss(s)).map(|(n, _)| n.as_str()).collect(); if !others.is_empty() { d.error( "model", @@ -521,7 +651,6 @@ impl Model { ), ); } - let axisymmetric = r.extent() == Some(Extent::Revolved); if axisymmetric && r.planar.values().any(|e| *e != Extent::Revolved) { d.error("model", "axisymmetric and plane (stress or strain) sections cannot be mixed in one model"); } @@ -542,6 +671,14 @@ impl Model { clockwise.push(label); } } + for sec in self.sections.values().filter(|s| truss(s)) { + for label in r.sets.get(&sec.set).map(|s| s.elements.as_slice()).unwrap_or_default() { + let Some(e) = self.elements.get(label) else { continue }; + if e.nodes.iter().filter_map(|n| self.nodes.get(n)).any(|n| n.xyz[2] != 0.0) { + off_plane.push(*label); + } + } + } if !off_plane.is_empty() { d.error("model", format!("2D elements {} have nodes off the XY plane (z ≠ 0); 2D models lie in the XY plane", first_labels(&off_plane))); } diff --git a/crates/model/src/constraints.rs b/crates/model/src/constraints.rs index 764d2b2..e0bda58 100644 --- a/crates/model/src/constraints.rs +++ b/crates/model/src/constraints.rs @@ -104,6 +104,7 @@ impl Model { .map(|(n, _)| format!("boundary condition `{n}`")), ); users.extend(self.loads.iter().filter(|(_, l)| l.reference_point() == Some(name)).map(|(n, _)| format!("load `{n}`"))); + users.extend(self.discrete_users(|t| t.reference_point() == Some(name))); not_in_use("reference point", name, users)?; remove_existing(&mut self.reference_points, "reference point", name).map(|_| ()) } diff --git a/crates/model/src/discrete.rs b/crates/model/src/discrete.rs new file mode 100644 index 0000000..a254b19 --- /dev/null +++ b/crates/model/src/discrete.rs @@ -0,0 +1,501 @@ +//! Discrete elements: point masses with rotary inertia, springs and dashpots. +//! +//! They are definitions, not mesh elements: each acts on a target (the nodes of a set, or a +//! reference point) and survives remeshing like a load does. Every node of the target gets its +//! own copy, with the full value: a point mass of 2 on a set of five nodes adds 2 at each of them +//! (a mass shared over a face is a point mass on a reference point coupled to the face). +//! +//! - A **point mass** adds a translational mass, the same in every direction or one per global +//! direction, and a rotary inertia tensor about the node (which needs rotational degrees of +//! freedom: a reference point, or a shell or beam node). +//! - A **spring** resists the relative displacement of its two ends along a direction with a +//! force `stiffness × relative displacement`; a **dashpot** resists the relative velocity with +//! `coefficient × relative velocity`. With no second end, the spring or dashpot connects each +//! node of the target to the ground. The direction is a global degree of freedom (1–3 +//! translations, 4–6 rotations: then a torsional spring or dashpot), a fixed vector, or the +//! line through the two ends (two-point only, from the undeformed positions). +//! +//! The linear (small displacement) form is the textbook spring element: for a unit direction n, +//! the element matrix is k [n nᵀ, −n nᵀ; −n nᵀ, n nᵀ] on the two ends' translations (R. D. Cook, +//! D. S. Malkus, M. E. Plesha, R. J. Witt, *Concepts and Applications of Finite Element +//! Analysis*, 4th ed., Wiley 2002, §2.2 and §2.10), and k n nᵀ for a grounded one. Dashpots +//! have the same form on velocities; static and frequency steps ignore them. + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::bcs::Target; +use crate::{Error, Model, Result, context, finite, insert_new, remove_existing, replace_existing}; + +/// Most copies of discrete elements a model resolves to (definitions × target nodes). +pub const MAX_DISCRETE: usize = 10_000_000; + +/// A translational mass: the same in every direction, or one per global direction. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(untagged)] +pub enum Mass { + Isotropic(f64), + Directional([f64; 3]), +} + +impl Mass { + /// The mass in X, Y and Z. + pub fn components(self) -> [f64; 3] { + match self { + Mass::Isotropic(m) => [m; 3], + Mass::Directional(m) => m, + } + } +} + +/// A point mass and rotary inertia at every node of a target. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct PointMass { + pub target: Target, + #[serde(default, skip_serializing_if = "Option::is_none")] + pub mass: Option, + /// The inertia tensor about the node in global axes, `[I11, I22, I33, I12, I13, I23]`: the + /// components of the symmetric tensor J with angular momentum H = J ω (products of inertia + /// carry their sign: I12 = −∫ x y dm). + #[serde(default, skip_serializing_if = "Option::is_none")] + pub inertia: Option<[f64; 6]>, +} + +impl PointMass { + pub fn validate(&self) -> Result<()> { + if self.mass.is_none() && self.inertia.is_none() { + return Err(Error::Invalid("a point mass needs a mass, an inertia or both".into())); + } + if let Some(m) = self.mass { + for c in m.components() { + if finite("mass", c)? < 0.0 { + return Err(Error::Invalid(format!("a mass is not negative, got {c}"))); + } + } + if m.components().iter().all(|c| *c == 0.0) && self.inertia.is_none() { + return Err(Error::Invalid("a point mass of zero adds nothing".into())); + } + } + if let Some(i) = self.inertia { + for c in i { + finite("inertia", c)?; + } + // Positive semidefinite: every principal minor is non-negative (to round-off). + let j = inertia_matrix(i); + let scale = i.iter().fold(0.0f64, |m, c| m.max(c.abs())); + let tolerance = 1e-12 * scale; + let minors = [ + j[0][0], + j[1][1], + j[2][2], + j[0][0] * j[1][1] - j[0][1] * j[0][1], + j[0][0] * j[2][2] - j[0][2] * j[0][2], + j[1][1] * j[2][2] - j[1][2] * j[1][2], + j[0][0] * (j[1][1] * j[2][2] - j[1][2] * j[2][1]) - j[0][1] * (j[1][0] * j[2][2] - j[1][2] * j[2][0]) + + j[0][2] * (j[1][0] * j[2][1] - j[1][1] * j[2][0]), + ]; + let powers = [1, 1, 1, 2, 2, 2, 3]; + if minors.iter().zip(powers).any(|(m, p)| *m < -tolerance * scale.powi(p - 1)) { + return Err(Error::Invalid( + "an inertia tensor must be positive semidefinite (non-negative moments, products no larger than the moments allow)".into(), + )); + } + if i.iter().all(|c| *c == 0.0) && self.mass.is_none_or(|m| m.components().iter().all(|c| *c == 0.0)) { + return Err(Error::Invalid("a point mass of zero adds nothing".into())); + } + } + Ok(()) + } +} + +/// `[I11, I22, I33, I12, I13, I23]` as a symmetric matrix. +pub fn inertia_matrix(i: [f64; 6]) -> [[f64; 3]; 3] { + let [i11, i22, i33, i12, i13, i23] = i; + [[i11, i12, i13], [i12, i22, i23], [i13, i23, i33]] +} + +/// The direction a spring or dashpot acts in. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum Direction { + /// A global degree of freedom: 1–3 translations along X, Y, Z, 4–6 rotations about them. + Dof(u8), + /// A fixed direction in global axes (translations; normalized when used). + Vector([f64; 3]), + /// The line through the two ends' undeformed positions (two-point only). + Line, +} + +impl Direction { + /// Whether it acts on rotations. + pub fn is_rotational(self) -> bool { + matches!(self, Direction::Dof(4..=6)) + } + + fn validate(self, two_point: bool) -> Result<()> { + match self { + Direction::Dof(d) if !(1..=6).contains(&d) => Err(Error::Invalid(format!("degree of freedom {d} is not in 1..=6"))), + Direction::Dof(_) => Ok(()), + Direction::Vector(v) => { + for c in v { + finite("direction", c)?; + } + if v.iter().all(|c| *c == 0.0) { + return Err(Error::Invalid("a direction vector must not be zero".into())); + } + Ok(()) + } + Direction::Line if !two_point => Err(Error::Invalid("the direction `line` needs a second end (it runs between the two ends)".into())), + Direction::Line => Ok(()), + } + } +} + +/// A linear spring or dashpot between the nodes of `first` and the ground, or between `first` +/// and `second` (each one node: a vertex set or a reference point). +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] +pub enum Spring { + /// Force (or moment) = stiffness × relative displacement (or rotation). + Spring { + first: Target, + #[serde(default, skip_serializing_if = "Option::is_none")] + second: Option, + direction: Direction, + stiffness: f64, + }, + /// Force (or moment) = coefficient × relative velocity (or angular velocity); dynamic steps + /// only. + Dashpot { + first: Target, + #[serde(default, skip_serializing_if = "Option::is_none")] + second: Option, + direction: Direction, + coefficient: f64, + }, +} + +impl Spring { + pub fn first(&self) -> &Target { + match self { + Spring::Spring { first, .. } | Spring::Dashpot { first, .. } => first, + } + } + + pub fn second(&self) -> Option<&Target> { + match self { + Spring::Spring { second, .. } | Spring::Dashpot { second, .. } => second.as_ref(), + } + } + + pub fn direction(&self) -> Direction { + match self { + Spring::Spring { direction, .. } | Spring::Dashpot { direction, .. } => *direction, + } + } + + /// The stiffness or the damping coefficient. + pub fn value(&self) -> f64 { + match self { + Spring::Spring { stiffness, .. } => *stiffness, + Spring::Dashpot { coefficient, .. } => *coefficient, + } + } + + pub fn is_dashpot(&self) -> bool { + matches!(self, Spring::Dashpot { .. }) + } + + /// `spring` or `dashpot`. + pub fn kind(&self) -> &'static str { + if self.is_dashpot() { "dashpot" } else { "spring" } + } + + /// The targets it acts on. + pub fn targets(&self) -> impl Iterator { + std::iter::once(self.first()).chain(self.second()) + } + + pub fn validate(&self) -> Result<()> { + let what = if self.is_dashpot() { "damping coefficient" } else { "stiffness" }; + let v = finite(what, self.value())?; + if v <= 0.0 { + return Err(Error::Invalid(format!("a {} needs a positive {what}, got {v}", self.kind()))); + } + self.direction().validate(self.second().is_some())?; + if self.second() == Some(self.first()) { + return Err(Error::Invalid(format!("a {} cannot connect a target to itself", self.kind()))); + } + Ok(()) + } +} + +/// A resolved direction: a global degree of freedom (1–6), or a unit vector acting on the +/// translations (from a fixed vector, or along the line of the ends). +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum ResolvedDirection { + Dof(u8), + Vector([f64; 3]), + Line([f64; 3]), +} + +impl ResolvedDirection { + /// The coefficients of the relative motion on each end's six degrees of freedom: the + /// spring's elongation is Σ c (u_second − u_first). + pub fn coefficients(self) -> [f64; 6] { + let mut c = [0.0; 6]; + match self { + ResolvedDirection::Dof(d) => { + if let Some(x) = usize::from(d).checked_sub(1).and_then(|i| c.get_mut(i)) { + *x = 1.0; + } + } + ResolvedDirection::Vector(n) | ResolvedDirection::Line(n) => c[..3].copy_from_slice(&n), + } + c + } + + pub fn is_rotational(self) -> bool { + matches!(self, ResolvedDirection::Dof(4..=6)) + } + + /// The unit vector it acts along (for a rotation, the axis it turns about). + pub fn axis(self) -> [f64; 3] { + let c = self.coefficients(); + if self.is_rotational() { [c[3], c[4], c[5]] } else { [c[0], c[1], c[2]] } + } +} + +/// One point mass at one node. +#[derive(Debug, Clone, PartialEq)] +pub struct ResolvedMass { + /// The definition it comes from. + pub name: String, + pub node: u32, + /// Mass in X, Y, Z. + pub mass: [f64; 3], + /// `[I11, I22, I33, I12, I13, I23]` (zero when there is none). + pub inertia: [f64; 6], +} + +impl ResolvedSpring { + /// `spring` or `dashpot`. + pub fn kind(&self) -> &'static str { + if self.dashpot { "dashpot" } else { "spring" } + } +} + +/// One spring or dashpot at one node (grounded) or between two nodes. +#[derive(Debug, Clone, PartialEq)] +pub struct ResolvedSpring { + pub name: String, + pub dashpot: bool, + pub first: u32, + /// `None`: to the ground. + pub second: Option, + pub direction: ResolvedDirection, + /// Stiffness, or damping coefficient. + pub value: f64, +} + +impl Model { + pub fn masses(&self) -> &BTreeMap { + &self.masses + } + + pub fn mass(&self, name: &str) -> Option<&PointMass> { + self.masses.get(name) + } + + pub fn springs(&self) -> &BTreeMap { + &self.springs + } + + pub fn spring(&self, name: &str) -> Option<&Spring> { + self.springs.get(name) + } + + fn check_mass(&self, name: &str, m: &PointMass) -> Result<()> { + m.validate().map_err(|e| context("point mass", name, e))?; + m.target.check_exists(self).map_err(|e| context("point mass", name, e)) + } + + pub fn add_mass(&mut self, name: &str, mass: PointMass) -> Result<()> { + self.check_mass(name, &mass)?; + insert_new(&mut self.masses, "point mass", name, mass) + } + + pub fn replace_mass(&mut self, name: &str, mass: PointMass) -> Result<()> { + self.check_mass(name, &mass)?; + replace_existing(&mut self.masses, "point mass", name, mass) + } + + pub fn remove_mass(&mut self, name: &str) -> Result<()> { + remove_existing(&mut self.masses, "point mass", name).map(|_| ()) + } + + fn check_spring(&self, name: &str, s: &Spring) -> Result<()> { + let what = s.kind(); + s.validate().map_err(|e| context(what, name, e))?; + for t in s.targets() { + t.check_exists(self).map_err(|e| context(what, name, e))?; + } + Ok(()) + } + + pub fn add_spring(&mut self, name: &str, spring: Spring) -> Result<()> { + self.check_spring(name, &spring)?; + insert_new(&mut self.springs, "spring", name, spring) + } + + pub fn replace_spring(&mut self, name: &str, spring: Spring) -> Result<()> { + self.check_spring(name, &spring)?; + replace_existing(&mut self.springs, "spring", name, spring) + } + + pub fn remove_spring(&mut self, name: &str) -> Result<()> { + remove_existing(&mut self.springs, "spring", name).map(|_| ()) + } + + /// The definitions among point masses, springs and dashpots that act on a target. + pub(crate) fn discrete_users(&self, uses: impl Fn(&Target) -> bool) -> Vec { + let mut users: Vec = self.masses.iter().filter(|(_, m)| uses(&m.target)).map(|(n, _)| format!("point mass `{n}`")).collect(); + users.extend(self.springs.iter().filter(|(_, s)| s.targets().any(&uses)).map(|(n, s)| format!("{} `{n}`", s.kind()))); + users + } + + /// The coordinates of a mesh node or a reference point's node label. + pub fn label_xyz(&self, label: u32) -> Option<[f64; 3]> { + self.nodes.get(&label).map(|n| n.xyz).or_else(|| self.reference_points.values().find(|p| p.label == label).map(|p| p.xyz)) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::Set; + + #[test] + fn definitions_validate_and_round_trip() { + let mut m = Model::default(); + m.add_set("tip", Set { nodes: vec![1], ..Set::default() }).unwrap(); + let json = serde_json::json!({"target": {"set": "tip"}, "mass": [1.0, 2.0, 0.0], "inertia": [1.0, 2.0, 3.0, 0.5, 0.0, 0.0]}); + let pm: PointMass = serde_json::from_value(json.clone()).unwrap(); + assert_eq!(pm.mass, Some(Mass::Directional([1.0, 2.0, 0.0]))); + assert_eq!(serde_json::to_value(&pm).unwrap(), json); + m.add_mass("pm", pm).unwrap(); + let iso: PointMass = serde_json::from_value(serde_json::json!({"target": {"set": "tip"}, "mass": 3})).unwrap(); + assert_eq!(iso.mass.map(Mass::components), Some([3.0; 3])); + for bad in [ + serde_json::json!({"target": {"set": "tip"}}), + serde_json::json!({"target": {"set": "tip"}, "mass": -1}), + serde_json::json!({"target": {"set": "tip"}, "mass": 0}), + serde_json::json!({"target": {"set": "tip"}, "inertia": [1, 1, 1, 2, 0, 0]}), + serde_json::json!({"target": {"set": "tip"}, "inertia": [-1, 1, 1, 0, 0, 0]}), + serde_json::json!({"target": {"set": "tip"}, "inertia": [0, 0, 0, 0, 0, 0]}), + ] { + let pm: PointMass = serde_json::from_value(bad.clone()).unwrap(); + assert!(m.add_mass("x", pm).is_err(), "{bad}"); + } + assert!(serde_json::from_value::(serde_json::json!({"target": {"set": "tip"}, "mass": 1, "masss": 2})).is_err()); + let missing: PointMass = serde_json::from_value(serde_json::json!({"target": {"set": "nope"}, "mass": 1})).unwrap(); + assert!(matches!(m.add_mass("x", missing), Err(Error::NotFound(_)))); + + let s: Spring = + serde_json::from_value(serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": {"dof": 3}, "stiffness": 100})) + .unwrap(); + m.add_spring("k", s).unwrap(); + let d: Spring = serde_json::from_value( + serde_json::json!({"type": "dashpot", "first": {"set": "tip"}, "direction": {"vector": [1, 1, 0]}, "coefficient": 2}), + ) + .unwrap(); + assert!(d.is_dashpot()); + m.add_spring("c", d).unwrap(); + for bad in [ + serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": "line", "stiffness": 1}), + serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": {"dof": 7}, "stiffness": 1}), + serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": {"vector": [0, 0, 0]}, "stiffness": 1}), + serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": {"dof": 1}, "stiffness": 0}), + serde_json::json!({"type": "spring", "first": {"set": "tip"}, "direction": {"dof": 1}, "stiffness": f64::NAN}), + serde_json::json!({"type": "dashpot", "first": {"set": "tip"}, "second": {"set": "tip"}, "direction": {"dof": 1}, "coefficient": 1}), + ] { + let Ok(s) = serde_json::from_value::(bad.clone()) else { continue }; + assert!(m.add_spring("x", s).is_err(), "{bad}"); + } + assert!( + serde_json::from_value::(serde_json::json!({"type": "dashpot", "first": {"set": "tip"}, "direction": "line", "stiffness": 1})) + .is_err() + ); + // The set is in use. + assert!(matches!(m.remove_set("tip"), Err(Error::InUse(_)))); + m.remove_mass("pm").unwrap(); + m.remove_spring("k").unwrap(); + m.remove_spring("c").unwrap(); + m.remove_set("tip").unwrap(); + } + + #[test] + fn directions_give_unit_coefficients() { + assert_eq!(ResolvedDirection::Dof(5).coefficients(), [0.0, 0.0, 0.0, 0.0, 1.0, 0.0]); + assert_eq!(ResolvedDirection::Line([0.6, 0.8, 0.0]).coefficients(), [0.6, 0.8, 0.0, 0.0, 0.0, 0.0]); + assert!(ResolvedDirection::Dof(4).is_rotational() && !ResolvedDirection::Dof(3).is_rotational()); + } + + #[test] + fn resolution_expands_and_checks_targets() { + use crate::{Procedure, ReferencePoint, Step}; + let mut m = Model::default(); + m.set_node(1, [0.0; 3]).unwrap(); + m.set_node(2, [3.0, 4.0, 0.0]).unwrap(); + m.add_set("a", Set { nodes: vec![1], ..Set::default() }).unwrap(); + m.add_set("b", Set { nodes: vec![2], ..Set::default() }).unwrap(); + m.add_set("ab", Set { nodes: vec![1, 2], ..Set::default() }).unwrap(); + m.add_reference_point("RP", ReferencePoint { xyz: [0.0, 0.0, 1.0], label: 0 }).unwrap(); + m.add_step(Step::new("Modes", Procedure::Frequency { modes: 1, min_hz: None, max_hz: None }), None).unwrap(); + let spring = |first: &str, second: Option<&str>, direction| Spring::Spring { + first: Target::Set(first.into()), + second: second.map(|s| Target::Set(s.into())), + direction, + stiffness: 10.0, + }; + m.add_spring("line", spring("a", Some("b"), Direction::Line)).unwrap(); + m.add_spring("ground", spring("ab", None, Direction::Vector([0.0, 0.0, 2.0]))).unwrap(); + m.add_mass("m", PointMass { target: Target::Set("ab".into()), mass: Some(Mass::Isotropic(2.0)), inertia: None }).unwrap(); + m.add_mass("j", PointMass { target: Target::ReferencePoint("RP".into()), mass: None, inertia: Some([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) }) + .unwrap(); + // No mesh elements: the discrete elements are the model. + let r = m.resolve().unwrap(); + assert_eq!(r.springs.len(), 3); + assert_eq!(r.springs[0].direction, ResolvedDirection::Vector([0.0, 0.0, 1.0])); + assert_eq!((r.springs[2].first, r.springs[2].second), (1, Some(2))); + assert_eq!(r.springs[2].direction, ResolvedDirection::Line([0.6, 0.8, 0.0])); + assert_eq!(r.masses.len(), 3); + assert_eq!(r.masses.iter().map(|x| x.node).collect::>(), vec![99_999_999, 1, 2]); + + // Errors: two-point ends with several nodes, rotations at solid nodes, inertia at mesh nodes. + m.add_spring("wide", spring("ab", Some("b"), Direction::Dof(1))).unwrap(); + m.add_spring("twist", spring("a", None, Direction::Dof(6))).unwrap(); + m.add_mass("bad", PointMass { target: Target::Set("a".into()), mass: None, inertia: Some([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) }).unwrap(); + let errors = m.resolve().unwrap_err().to_string(); + for want in ["spring `wide`: a two-point spring connects one node to one node", "spring `twist`: acts on a rotation", "point mass `bad`"] { + assert!(errors.contains(want), "{errors}"); + } + // A reference point used by a discrete element cannot be deleted. + assert!(matches!(m.remove_reference_point("RP"), Err(Error::InUse(_)))); + // Torsion at a reference point is fine. + m.remove_spring("twist").unwrap(); + m.remove_spring("wide").unwrap(); + m.remove_mass("bad").unwrap(); + m.add_spring( + "twist", + Spring::Spring { first: Target::ReferencePoint("RP".into()), second: None, direction: Direction::Dof(6), stiffness: 1.0 }, + ) + .unwrap(); + m.resolve().unwrap(); + // The definitions survive a file round trip. + let back: Model = serde_json::from_str(&serde_json::to_string(&m).unwrap()).unwrap(); + assert_eq!(back, m); + } +} diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index fa64148..5dec6cc 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -8,7 +8,8 @@ //! - named sets ([`sets`]) defined by geometry names or explicit labels, resolved to nodes, //! elements and element faces against the current mesh, so they survive remeshing; //! - materials and the unit system ([`materials`], [`units`]), sections with their element -//! formulation ([`sections`]), reference points and couplings ([`constraints`]); +//! formulation ([`sections`]), reference points and couplings ([`constraints`]), point masses, +//! springs and dashpots ([`discrete`]); //! - boundary conditions ([`bcs`]), loads ([`loads`]) and steps with an explicit table of which //! boundary conditions and loads are active in each step ([`steps`]); //! - the check that names every missing or bad item, and the resolved view the deck writers and @@ -33,6 +34,7 @@ pub mod bcs; pub mod check; pub mod constraints; +pub mod discrete; pub mod fluids; pub mod geometry; pub mod loads; @@ -51,6 +53,7 @@ use serde::{Deserialize, Serialize}; pub use bcs::{BoundaryCondition, Dofs, Symmetry, Target}; pub use check::{Diagnostic, Resolved, ResolvedSet, Severity}; pub use constraints::{Coupling, CouplingKind, ReferencePoint}; +pub use discrete::{Direction, Mass, PointMass, ResolvedDirection, ResolvedMass, ResolvedSpring, Spring}; pub use fluids::{BoxFace, BoxFaces, FaceCondition, FluidBoundary, FluidDomain, Profile, SurfaceBoundary, SurfaceCondition}; pub use geometry::{Classification, Entity, GeometryRef}; pub use loads::{Distribution, Load, centrifugal_acceleration}; @@ -113,6 +116,10 @@ pub struct Model { sections: BTreeMap, reference_points: BTreeMap, couplings: BTreeMap, + #[serde(skip_serializing_if = "BTreeMap::is_empty")] + masses: BTreeMap, + #[serde(skip_serializing_if = "BTreeMap::is_empty")] + springs: BTreeMap, boundary_conditions: BTreeMap, loads: BTreeMap, #[serde(skip_serializing_if = "BTreeMap::is_empty")] @@ -144,6 +151,10 @@ struct ModelData { #[serde(default)] couplings: BTreeMap, #[serde(default)] + masses: BTreeMap, + #[serde(default)] + springs: BTreeMap, + #[serde(default)] boundary_conditions: BTreeMap, #[serde(default)] loads: BTreeMap, @@ -198,6 +209,12 @@ impl TryFrom for Model { for (n, v) in d.couplings { m.add_coupling(&n, v)?; } + for (n, v) in d.masses { + m.add_mass(&n, v)?; + } + for (n, v) in d.springs { + m.add_spring(&n, v)?; + } for (n, v) in d.boundary_conditions { m.add_boundary_condition(&n, v)?; } @@ -307,6 +324,8 @@ impl Model { sections: self.sections.clone(), reference_points: self.reference_points.clone(), couplings: self.couplings.clone(), + masses: self.masses.clone(), + springs: self.springs.clone(), boundary_conditions: self.boundary_conditions.clone(), loads: self.loads.clone(), fluid_domains: self.fluid_domains.clone(), @@ -315,13 +334,15 @@ impl Model { } /// Counts of each kind of definition, for summaries. - pub fn definition_counts(&self) -> [(&'static str, usize); 9] { + pub fn definition_counts(&self) -> [(&'static str, usize); 11] { [ ("sets", self.sets.len()), ("materials", self.materials.len()), ("sections", self.sections.len()), ("reference_points", self.reference_points.len()), ("couplings", self.couplings.len()), + ("masses", self.masses.len()), + ("springs", self.springs.len()), ("boundary_conditions", self.boundary_conditions.len()), ("loads", self.loads.len()), ("fluid_domains", self.fluid_domains.len()), diff --git a/crates/model/src/mesh.rs b/crates/model/src/mesh.rs index 019ae13..c1552cc 100644 --- a/crates/model/src/mesh.rs +++ b/crates/model/src/mesh.rs @@ -14,6 +14,7 @@ pub struct Node { /// Element topology. Node order is corners first, then mid-side nodes, with positive volume for /// a right-handed ordering: +/// - rod3: the two ends, then the middle node; /// - tet10: mid-side nodes 5–10 on edges 1-2, 2-3, 3-1, 1-4, 2-4, 3-4; /// - wedge6: triangle 1-2-3, then 4-5-6 on the side its right-hand normal points to; /// - hex8/hex20: face 1-2-3-4, then 5-6-7-8 above it; hex20's mid-side nodes 9–12 lie on the @@ -25,8 +26,10 @@ pub struct Node { #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] #[serde(rename_all = "snake_case")] pub enum ElementKind { - /// Two-node axial member. + /// Two-node axial member (truss). Rod2, + /// Three-node axial member (quadratic truss, straight or curved). + Rod3, /// Two-node beam. Beam2, /// Three-node triangle (shell or membrane). @@ -50,8 +53,9 @@ pub enum ElementKind { } impl ElementKind { - pub const ALL: [ElementKind; 11] = [ + pub const ALL: [ElementKind; 12] = [ ElementKind::Rod2, + ElementKind::Rod3, ElementKind::Beam2, ElementKind::Tri3, ElementKind::Tri6, @@ -67,7 +71,7 @@ impl ElementKind { pub fn node_count(self) -> usize { match self { ElementKind::Rod2 | ElementKind::Beam2 => 2, - ElementKind::Tri3 => 3, + ElementKind::Rod3 | ElementKind::Tri3 => 3, ElementKind::Quad4 | ElementKind::Tet4 => 4, ElementKind::Tri6 | ElementKind::Wedge6 => 6, ElementKind::Quad8 | ElementKind::Hex8 => 8, @@ -79,7 +83,7 @@ impl ElementKind { /// Corner nodes: the first ones in the node list (the rest are mid-side nodes). pub fn corner_count(self) -> usize { match self { - ElementKind::Rod2 | ElementKind::Beam2 => 2, + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => 2, ElementKind::Tri3 | ElementKind::Tri6 => 3, ElementKind::Quad4 | ElementKind::Quad8 | ElementKind::Tet4 | ElementKind::Tet10 => 4, ElementKind::Wedge6 => 6, @@ -95,7 +99,7 @@ impl ElementKind { /// Topological dimension: 1 for line elements, 2 for shells, 3 for solids. pub fn dimension(self) -> u8 { match self { - ElementKind::Rod2 | ElementKind::Beam2 => 1, + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => 1, ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => 2, ElementKind::Tet4 | ElementKind::Tet10 | ElementKind::Wedge6 | ElementKind::Hex8 | ElementKind::Hex20 => 3, } @@ -105,6 +109,7 @@ impl ElementKind { pub fn name(self) -> &'static str { match self { ElementKind::Rod2 => "rod2", + ElementKind::Rod3 => "rod3", ElementKind::Beam2 => "beam2", ElementKind::Tri3 => "tri3", ElementKind::Tri6 => "tri6", @@ -134,7 +139,7 @@ impl ElementKind { /// - Line elements have none. pub fn faces(self) -> &'static [&'static [usize]] { match self { - ElementKind::Rod2 | ElementKind::Beam2 => &[], + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => &[], ElementKind::Tri3 => &[&[0, 1, 2], &[0, 2, 1]], ElementKind::Tri6 => &[&[0, 1, 2, 3, 4, 5], &[0, 2, 1, 5, 4, 3]], ElementKind::Quad4 => &[&[0, 1, 2, 3], &[0, 3, 2, 1]], @@ -160,6 +165,7 @@ impl ElementKind { pub fn edges(self) -> &'static [&'static [usize]] { match self { ElementKind::Rod2 | ElementKind::Beam2 => &[&[0, 1]], + ElementKind::Rod3 => &[&[0, 1, 2]], ElementKind::Tri3 => &[&[0, 1], &[1, 2], &[2, 0]], ElementKind::Tri6 => &[&[0, 1, 3], &[1, 2, 4], &[2, 0, 5]], ElementKind::Quad4 => &[&[0, 1], &[1, 2], &[2, 3], &[3, 0]], @@ -357,7 +363,7 @@ mod tests { /// Reference coordinates with positive volume for each kind (mid-side nodes at edge midpoints). fn reference(kind: ElementKind) -> Vec<[f64; 3]> { let corners: Vec<[f64; 3]> = match kind { - ElementKind::Rod2 | ElementKind::Beam2 => vec![[0., 0., 0.], [1., 0., 0.]], + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => vec![[0., 0., 0.], [1., 0., 0.]], ElementKind::Tri3 | ElementKind::Tri6 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]], ElementKind::Quad4 | ElementKind::Quad8 => vec![[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.]], ElementKind::Tet4 | ElementKind::Tet10 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]], diff --git a/crates/model/src/sets.rs b/crates/model/src/sets.rs index 4ee451f..2df63c7 100644 --- a/crates/model/src/sets.rs +++ b/crates/model/src/sets.rs @@ -122,6 +122,7 @@ impl Model { let mut users = Vec::new(); users.extend(self.sections.iter().filter(|(_, s)| s.set == set).map(|(n, _)| format!("section `{n}`"))); users.extend(self.couplings.iter().filter(|(_, c)| c.set == set).map(|(n, _)| format!("coupling `{n}`"))); + users.extend(self.discrete_users(|t| t.set() == Some(set))); users .extend(self.boundary_conditions.iter().filter(|(_, b)| b.target().set() == Some(set)).map(|(n, _)| format!("boundary condition `{n}`"))); users.extend(self.loads.iter().filter(|(_, l)| l.set() == Some(set)).map(|(n, _)| format!("load `{n}`"))); diff --git a/crates/render/src/display.rs b/crates/render/src/display.rs index a66e294..d3630d0 100644 --- a/crates/render/src/display.rs +++ b/crates/render/src/display.rs @@ -79,7 +79,7 @@ pub struct DisplayMesh { /// counter-clockwise seen from outside for a positively oriented element. fn faces(kind: ElementKind) -> &'static [&'static [usize]] { match kind { - ElementKind::Rod2 | ElementKind::Beam2 => &[], + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => &[], ElementKind::Tri3 | ElementKind::Tri6 => &[&[0, 1, 2]], ElementKind::Quad4 | ElementKind::Quad8 => &[&[0, 1, 2, 3]], ElementKind::Tet4 | ElementKind::Tet10 => &[&[0, 2, 1], &[0, 1, 3], &[1, 2, 3], &[2, 0, 3]], diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 5f20529..91e7f46 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -12,7 +12,7 @@ use crate::constraints::{Budget, Reduction, coupling_rows}; use crate::eigen::{Operator, Wanted, lanczos}; use crate::ordering::nested_dissection; use crate::sparse::{Csr, NodeGraph, add_block, add_diagonal, add_scalar}; -use crate::system::{Mesh, inelastic_warnings, loads, prescribed}; +use crate::system::{Mesh, inelastic_warnings, loads, prescribed, translation_block}; use crate::{Error, Limits, Result}; /// Size and effort of a solution. @@ -147,6 +147,24 @@ fn assemble(mesh: &Mesh, graph: &NodeGraph, mass: Option) -> Result add_scalar(mm, layout, &t.nodes, kinds, &truss.mass(rho, t.area).map_err(at)?.data)?, + MassKind::Lumped => add_diagonal(mm, layout, &t.nodes, kinds, &truss.lumped_mass(rho, t.area).map_err(at)?)?, + } + } + } + mesh.discrete.add_stiffness(&mut k)?; + if let Some(mm) = m.as_mut() { + mesh.discrete.add_mass(mm)?; + } Ok(Assembled { k, m }) } @@ -184,6 +202,32 @@ fn split(full: &[f64], mesh: &Mesh) -> (Vec<[f64; 3]>, Vec<[f64; 6]>) { ((0..l.nodes()).map(|n| directions(full, l, n)).collect(), (0..mesh.points.len()).map(|k| directions(full, l, l.point_slot(k))).collect()) } +/// The matrix pattern of the mesh (node blocks of its elements, trusses and discrete elements) +/// and the order of its nodes for the factorization. +fn pattern(mesh: &Mesh) -> (NodeGraph, Vec) { + let graph = NodeGraph::new(mesh.block_count(), mesh.blocks()); + // The order covers the mesh nodes; reference point equations follow them. + let n = mesh.labels.len(); + let order = if graph.adjacency.len() == n { + nested_dissection(&mesh.xyz, &graph.adjacency) + } else { + let mesh_part: Vec> = graph.adjacency.iter().take(n).map(|row| row.iter().copied().filter(|b| *b < n).collect()).collect(); + nested_dissection(&mesh.xyz, &mesh_part) + }; + (graph, order) +} + +/// Warnings for what a step without velocities ignores. +fn ignored_dashpots(mesh: &Mesh, warnings: &mut Vec) { + let names = mesh.discrete.dashpots(); + if !names.is_empty() { + warnings.push(format!( + "dashpots {} act in dynamic steps only; this step ignores them", + names.iter().map(|n| format!("`{n}`")).collect::>().join(", ") + )); + } +} + /// Constraint rows, prescribed values and their elimination for a step. /// The model, its resolved view, the mesh and the step being solved. struct Setup<'a> { @@ -210,8 +254,8 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).collect(); warnings.extend(inelastic_warnings(model)); - let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice())); - let order = nested_dissection(&mesh.xyz, &graph.adjacency); + let (graph, order) = pattern(&mesh); + ignored_dashpots(&mesh, &mut warnings); let setup = Setup { model, resolved: &resolved, mesh: &mesh, step }; let red = reduction(&setup, true, &order, budget(limits, &mesh, &graph, 1), &mut warnings)?; let f = loads(model, &resolved, &mesh, step, &red.active)?; @@ -276,10 +320,28 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul } elements.push(ElementStress { label: e.label, kind: e.kind, points, nodal }); } + for t in &mesh.trusses { + let truss = t.truss()?; + let at = |err| Error::element(t.label, err); + let ue: Vec<[f64; 3]> = t.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect(); + let points = truss.stresses(t.young, &ue).map_err(at)?; + let nodal = truss.extrapolate(&points.iter().map(|p| p.stress).collect::>()).map_err(at)?; + for (n, s) in t.nodes.iter().zip(&nodal) { + if let (Some(acc), Some(c)) = (sum.get_mut(*n), count.get_mut(*n)) { + for k in 0..6 { + acc[k] += s[k]; + } + *c += 1; + } + } + elements.push(ElementStress { label: t.label, kind: t.kind, points, nodal }); + } + elements.sort_by_key(|e| e.label); let stress: Vec<[f64; 6]> = sum.iter().zip(&count).map(|(s, c)| s.map(|v| if *c > 0 { v / f64::from(*c) } else { 0.0 })).collect(); let von_mises = stress.iter().map(von_mises).collect(); let scale = applied_total.iter().chain(&reaction_total).fold(0.0f64, |m, v| m.max(v.abs())); - let imbalance = (0..3).map(|c| (applied_total[c] + reaction_total[c]).abs()).fold(0.0, f64::max); + let ground = mesh.discrete.ground_force(&u, &mesh); + let imbalance = (0..3).map(|c| (applied_total[c] + reaction_total[c] + ground[c]).abs()).fold(0.0, f64::max); if scale > 0.0 && imbalance > 1e-6 * scale { warnings.push(format!("loads and reactions differ by {imbalance:.3e} (relative {:.1e})", imbalance / scale)); } @@ -338,14 +400,23 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).collect(); warnings.extend(inelastic_warnings(model)); - let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice())); - let order = nested_dissection(&mesh.xyz, &graph.adjacency); + let (graph, order) = pattern(&mesh); + ignored_dashpots(&mesh, &mut warnings); let setup = Setup { model, resolved: &resolved, mesh: &mesh, step }; let red = reduction(&setup, false, &order, budget(limits, &mesh, &graph, 2), &mut warnings)?; let Assembled { k, m } = assemble(&mesh, &graph, Some(options.mass))?; let m = m.ok_or_else(|| Error::Internal("mass matrix".into()))?; // Every element's mass entries sum to its mass, once per direction. - let total_mass = m.val.iter().sum::() / if mesh.theory.is_some() { 2.0 } else { 3.0 }; + // (the mean over the directions for a point mass given per direction; rotary inertia is left + // out), added in storage order. + let translation = |d: usize| mesh.kind(d).and_then(DofKind::direction).is_some_and(|c| c < 3); + let mut translational = 0.0; + for d in (0..m.n).filter(|d| translation(*d)) { + for v in m.row(d).1 { + translational += v; + } + } + let total_mass = translational / mesh.translation_kinds().len() as f64; let neq = red.neq(); if neq == 0 { return Err(Error::Invalid("every degree of freedom is held; there is nothing to vibrate".into())); diff --git a/crates/solver/src/constraints.rs b/crates/solver/src/constraints.rs index fd5ee5a..202a295 100644 --- a/crates/solver/src/constraints.rs +++ b/crates/solver/src/constraints.rs @@ -191,10 +191,15 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) -> } /// The degrees of freedom that exist: every mesh node's, and the reference point directions -/// some coupling row involves. +/// some coupling row or discrete element involves. fn active_dofs(mesh: &Mesh, rows: &[Row]) -> Vec { let mut active = vec![false; mesh.dofs()]; active.iter_mut().take(mesh.mesh_dofs()).for_each(|a| *a = true); + for dof in mesh.discrete.dofs() { + if let Some(a) = active.get_mut(dof) { + *a = true; + } + } for row in rows { for (dof, c) in &row.terms { if *c != 0.0 diff --git a/crates/solver/src/discrete.rs b/crates/solver/src/discrete.rs new file mode 100644 index 0000000..fe21b42 --- /dev/null +++ b/crates/solver/src/discrete.rs @@ -0,0 +1,198 @@ +//! Point masses, rotary inertia, springs and dashpots in the implicit system. +//! +//! The model resolves each definition to its nodes ([`modelcraft_model::Resolved::masses`] and +//! `springs`); here each copy is turned into global degrees of freedom and its matrix: +//! +//! - a point mass adds m_x, m_y, m_z on the node's translations and its inertia tensor J (3 × 3, +//! off-diagonal products included) on its rotations, the same in the consistent and the lumped +//! mass matrix; +//! - a spring with stiffness k and unit direction coefficients c on each end's six degrees of +//! freedom adds k c cᵀ to each end and −k c cᵀ between them (Cook, Malkus, Plesha, Witt 2002, +//! §2.2): the elongation is cᵀ(u_second − u_first); +//! - a dashpot has the same form on velocities; static and frequency steps have no velocities, so +//! they ignore dashpots with a warning (damping enters dynamic steps). +//! +//! What each copy needs at its nodes ([`needs`]) goes into the mesh's [`modelcraft_fem::DofLayout`] +//! like an element's declared kinds; [`dof_of`] then finds its unknowns there. + +use modelcraft_fem::{DofKind, DofSet}; +use modelcraft_model::{Resolved, ResolvedSpring}; + +use crate::sparse::{Csr, add_entry}; +use crate::system::Mesh; +use crate::{Error, Result}; + +/// A point mass on global degrees of freedom. +#[derive(Debug, Clone, PartialEq)] +pub struct MassAt { + pub name: String, + /// Translations along X, Y, Z: (degree of freedom, mass). + pub translations: Vec<(usize, f64)>, + /// The rotations' degrees of freedom and the inertia tensor on them (empty without inertia). + pub rotations: Vec, + pub inertia: [[f64; 3]; 3], +} + +/// A spring or dashpot on global degrees of freedom: the elongation is Σ c u over `terms`, whose +/// coefficients are the direction's on the second end and their negatives on the first. +#[derive(Debug, Clone, PartialEq)] +pub struct SpringAt { + pub name: String, + pub dashpot: bool, + /// Connected to the ground (one end only). + pub grounded: bool, + pub terms: Vec<(usize, f64)>, + pub value: f64, +} + +/// Every discrete element of the model on the mesh's degrees of freedom. +#[derive(Debug, Clone, Default, PartialEq)] +pub struct Discrete { + pub masses: Vec, + pub springs: Vec, +} + +/// The kinds each discrete copy needs at the labels it acts on: translations for a mass, +/// rotations for an inertia, translations or rotations for a spring or dashpot. +/// `translations` are the model's (X and Y in a 2D model). +pub fn needs(resolved: &Resolved, translations: DofSet) -> impl Iterator + '_ { + let masses = resolved.masses.iter().flat_map(move |m| { + let mass = m.mass.iter().any(|c| *c != 0.0).then_some((m.node, translations)); + let inertia = m.inertia.iter().any(|c| *c != 0.0).then_some((m.node, DofSet::ROTATIONS)); + mass.into_iter().chain(inertia) + }); + // A node of a 3D model moves in all three directions (a spring along X does not take Y and + // Z away), so a spring asks for every translation, or every rotation. + let springs = resolved.springs.iter().flat_map(move |s| { + let kinds = if s.direction.is_rotational() { DofSet::ROTATIONS } else { translations }; + std::iter::once((s.first, kinds)).chain(s.second.map(|b| (b, kinds))) + }); + masses.chain(springs) +} + +/// The global degree of freedom of direction `d` (0–2 translations, 3–5 rotations) of a node or +/// reference point label, or `None` when the node has no such direction. +pub fn dof_of(mesh: &Mesh, label: u32, d: usize) -> Option { + mesh.layout.dof(mesh.layout.slot(label)?, DofKind::from_direction(d)?) +} + +fn missing(what: &str, name: &str, label: u32, d: usize) -> Error { + let direction = DofKind::from_direction(d).map_or("?", DofKind::describe); + Error::Invalid(format!("{what} `{name}`: node {label} has no degree of freedom in {direction}")) +} + +/// The spring's terms at one end, with the direction's coefficients times `sign`. +fn ends(mesh: &Mesh, s: &ResolvedSpring, label: u32, sign: f64) -> Result> { + s.direction + .coefficients() + .iter() + .enumerate() + .filter(|(_, c)| **c != 0.0) + .map(|(d, c)| dof_of(mesh, label, d).map(|dof| (dof, sign * c)).ok_or_else(|| missing(s.kind(), &s.name, label, d))) + .collect() +} + +impl Discrete { + /// The resolved discrete elements on the mesh's degrees of freedom. + pub fn new(mesh: &Mesh, resolved: &Resolved) -> Result { + let mut out = Discrete::default(); + for m in &resolved.masses { + // The model's translations only (a 2D model's masses act in its plane). + let planar = mesh.translation_kinds(); + let translations: Vec<(usize, f64)> = (0..3) + .filter(|d| m.mass[*d] != 0.0 && DofKind::translation(*d).is_some_and(|k| planar.contains(k))) + .map(|d| dof_of(mesh, m.node, d).map(|dof| (dof, m.mass[d])).ok_or_else(|| missing("point mass", &m.name, m.node, d))) + .collect::>()?; + let inertia = modelcraft_model::discrete::inertia_matrix(m.inertia); + let rotations = if m.inertia.iter().any(|c| *c != 0.0) { + (3..6).map(|d| dof_of(mesh, m.node, d).ok_or_else(|| missing("point mass", &m.name, m.node, d))).collect::>()? + } else { + Vec::new() + }; + out.masses.push(MassAt { name: m.name.clone(), translations, rotations, inertia }); + } + for s in &resolved.springs { + let mut terms = ends(mesh, s, s.first, -1.0)?; + if let Some(b) = s.second { + terms.extend(ends(mesh, s, b, 1.0)?); + } + out.springs.push(SpringAt { name: s.name.clone(), dashpot: s.dashpot, grounded: s.second.is_none(), terms, value: s.value }); + } + Ok(out) + } + + /// The groups of degrees of freedom each element couples, for the matrix pattern. + pub fn couplings(&self) -> impl Iterator> + '_ { + let masses = self.masses.iter().map(|m| m.translations.iter().map(|t| t.0).chain(m.rotations.iter().copied()).collect()); + masses.chain(self.springs.iter().map(|s| s.terms.iter().map(|t| t.0).collect())) + } + + /// Every degree of freedom a discrete element acts on. + pub fn dofs(&self) -> impl Iterator + '_ { + self.couplings().flatten() + } + + /// Names of the dashpots, which static and frequency steps ignore. + pub fn dashpots(&self) -> Vec<&str> { + let mut names: Vec<&str> = self.springs.iter().filter(|s| s.dashpot).map(|s| s.name.as_str()).collect(); + names.dedup(); + names + } + + /// Add the springs' stiffness. + pub fn add_stiffness(&self, k: &mut Csr) -> Result<()> { + for s in self.springs.iter().filter(|s| !s.dashpot) { + for (i, ci) in &s.terms { + for (j, cj) in &s.terms { + add_entry(k, *i, *j, s.value * ci * cj)?; + } + } + } + Ok(()) + } + + /// Add the point masses and rotary inertia. + pub fn add_mass(&self, m: &mut Csr) -> Result<()> { + for p in &self.masses { + for (d, v) in &p.translations { + add_entry(m, *d, *d, *v)?; + } + for (a, da) in p.rotations.iter().enumerate() { + for (b, db) in p.rotations.iter().enumerate() { + let v = p.inertia[a][b]; + if v != 0.0 { + add_entry(m, *da, *db, v)?; + } + } + } + } + Ok(()) + } + + /// The sum of the forces springs to the ground exert on the model, along X, Y and Z (the + /// loads they carry to the ground, besides the boundary conditions' reactions). + pub fn ground_force(&self, u: &[f64], mesh: &Mesh) -> [f64; 3] { + let mut total = [0.0; 3]; + for s in self.springs.iter().filter(|s| !s.dashpot && s.grounded) { + // The spring's force on the model: −k (cᵀu) c over its terms. + let stretch: f64 = s.terms.iter().map(|(d, c)| c * u.get(*d).copied().unwrap_or(0.0)).sum(); + for (d, c) in &s.terms { + if let Some(t) = mesh.kind(*d).and_then(DofKind::direction).filter(|d| *d < 3).and_then(|d| total.get_mut(d)) { + *t -= s.value * stretch * c; + } + } + } + total + } + + /// Gravity on the point masses: m_d g_d on each translation. + pub fn add_gravity(&self, f: &mut [f64], g: [f64; 3], mesh: &Mesh) { + for p in &self.masses { + for (dof, m) in &p.translations { + if let (Some(x), Some(a)) = (f.get_mut(*dof), mesh.kind(*dof).and_then(DofKind::direction).and_then(|d| g.get(d))) { + *x += m * a; + } + } + } + } +} diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index 56fffc0..50db5b9 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -25,6 +25,7 @@ pub mod analysis; pub mod cholesky; pub mod constraints; +pub mod discrete; pub mod eigen; pub mod ordering; pub mod sparse; diff --git a/crates/solver/src/sparse.rs b/crates/solver/src/sparse.rs index c078d32..e57b5a3 100644 --- a/crates/solver/src/sparse.rs +++ b/crates/solver/src/sparse.rs @@ -103,10 +103,10 @@ impl NodeGraph { self.adjacency.iter().enumerate().map(|(a, row)| layout.range(a).len() * row.iter().map(|b| layout.range(*b).len()).sum::()).sum() } - /// A zero matrix over the mesh nodes' degrees of freedom (layout slots `0..nodes`) on this - /// pattern. + /// A zero matrix over the degrees of freedom of the graph's layout slots (the mesh nodes + /// `0..nodes`, and reference points when the graph has rows for them) on this pattern. pub fn zero_matrix(&self, layout: &DofLayout) -> Csr { - let n = layout.node_dofs(); + let n = self.adjacency.len().checked_sub(1).map_or(0, |last| layout.range(last).end); let mut ptr = Vec::with_capacity(n + 1); ptr.push(0); let mut idx = Vec::with_capacity(self.entries(layout)); @@ -210,6 +210,14 @@ pub fn add_diagonal(k: &mut Csr, layout: &DofLayout, nodes: &[usize], kinds: Dof Ok(()) } +/// Add `v` to entry (`i`, `j`), which must be in the pattern. +pub fn add_entry(k: &mut Csr, i: usize, j: usize, v: f64) -> Result<()> { + let (start, end) = (k.ptr.get(i).copied().unwrap_or(0), k.ptr.get(i + 1).copied().unwrap_or(0)); + let at = k.idx.get(start..end).and_then(|cols| cols.binary_search(&j).ok()).ok_or_else(|| Error::Internal("entry outside the pattern".into()))?; + add_at(k, start + at, v); + Ok(()) +} + #[cfg(test)] mod tests { use super::*; diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 1179394..943eaab 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -1,5 +1,6 @@ -//! From the resolved model to the discrete system: the solid elements with their materials and -//! formulations, the coupled reference points, the degrees of freedom, prescribed values and the +//! From the resolved model to the discrete system: the solid, 2D and truss elements with their +//! materials and formulations, point masses, springs and dashpots ([`crate::discrete`]), the +//! reference points couplings or discrete elements use, the degrees of freedom, prescribed values and the //! load vector of a step. //! //! Degrees of freedom ([`DofLayout`]): the kinds each mesh node's elements declare (three @@ -8,6 +9,7 @@ use modelcraft_fem::element::{Formulation, Solid}; use modelcraft_fem::plane::{Plane, Theory}; +use modelcraft_fem::truss::{Truss, young_modulus}; use modelcraft_fem::{DofKind, DofLayout, DofSet, Elasticity, ElementDofs, PointStress, elasticity, face, plane}; use modelcraft_model::nodal::edge_pressure; use std::collections::BTreeSet; @@ -117,6 +119,42 @@ impl Elem { } } +/// A truss element (rod2, rod3) ready for the kernel. +#[derive(Debug, Clone)] +pub struct Bar { + pub label: u32, + pub kind: ElementKind, + /// Indices into [`Mesh::labels`] (layout slots). + pub nodes: Vec, + pub xyz: Vec<[f64; 3]>, + pub young: f64, + pub area: f64, + pub density: Option, +} + +impl Bar { + pub fn truss(&self) -> Result> { + Truss::new(self.kind, &self.xyz).map_err(|e| Error::element(self.label, e)) + } +} + +/// The translations of a model's nodes: X and Y in a 2D model, else X, Y and Z. Trusses and +/// translational discrete elements declare these. +pub fn translation_kinds(planar: bool) -> DofSet { + if planar { DofSet::PLANE } else { DofSet::TRANSLATIONS } +} + +/// Keep the rows and columns of a matrix over three translations per node (node by node) that +/// `kinds` has: the in-plane block of a truss in a 2D model, or the matrix itself. +pub fn translation_block(data: &[f64], nodes: usize, kinds: DofSet) -> Vec { + if kinds == DofSet::TRANSLATIONS { + return data.to_vec(); + } + let keep: Vec = (0..3 * nodes).filter(|i| DofKind::translation(i % 3).is_some_and(|k| kinds.contains(k))).collect(); + let width = 3 * nodes; + keep.iter().flat_map(|i| keep.iter().map(move |j| data.get(i * width + j).copied().unwrap_or(0.0))).collect() +} + /// A reference point some coupling ties to the mesh. #[derive(Debug, Clone, PartialEq)] pub struct Point { @@ -138,6 +176,9 @@ pub struct Mesh { pub layout: DofLayout, /// The 2D theory when this is a 2D model (a model is 2D as a whole). pub theory: Option, + pub trusses: Vec, + /// Point masses, springs and dashpots ([`crate::discrete`]). + pub discrete: crate::discrete::Discrete, } impl Mesh { @@ -145,6 +186,15 @@ impl Mesh { self.labels.binary_search(&label).ok() } + pub fn truss(&self, label: u32) -> Option<&Bar> { + self.trusses.binary_search_by_key(&label, |e| e.label).ok().and_then(|i| self.trusses.get(i)) + } + + /// The kinds trusses and translational discrete elements carry at their nodes. + pub fn translation_kinds(&self) -> DofSet { + translation_kinds(self.theory.is_some()) + } + pub fn element(&self, label: u32) -> Option<&Elem> { self.elements.binary_search_by_key(&label, |e| e.label).ok().and_then(|i| self.elements.get(i)) } @@ -188,17 +238,32 @@ impl Mesh { /// Build it from the model and its resolved view, refusing what the solver cannot represent. pub fn new(model: &Model, resolved: &Resolved, max_dofs: usize) -> Result { let mut elements = Vec::new(); + let mut trusses = Vec::new(); // A model is 2D as a whole (the model check refuses a mix). let planar = model.sections().values().find_map(|s| Theory::from_section(&s.kind)); + if matches!(planar, Some(Theory::Axisymmetric)) + && let Some(what) = model + .sections() + .values() + .any(|s| matches!(s.kind, SectionKind::Rod { .. })) + .then_some("trusses") + .or_else(|| (!model.masses().is_empty() || !model.springs().is_empty()).then_some("point masses, springs and dashpots")) + { + return Err(Error::NotSupported(format!("{what} in axisymmetric models"))); + } for (name, section) in model.sections() { let theory = Theory::from_section(§ion.kind); - if !matches!(section.kind, SectionKind::Solid {}) && theory.is_none() { + let rod = match section.kind { + SectionKind::Rod { area } => Some(area), + _ => None, + }; + if !matches!(section.kind, SectionKind::Solid {}) && theory.is_none() && rod.is_none() { return Err(Error::NotSupported(format!( - "{} sections (section `{name}`); the solver has solid and 2D (plane stress, plane strain, axisymmetric) elements so far", + "{} sections (section `{name}`); the solver has solid, 2D (plane stress, plane strain, axisymmetric) and truss elements so far", section.kind.name() ))); } - if theory.is_some() != planar.is_some() { + if theory.is_some() != planar.is_some() && rod.is_none() { return Err(Error::NotSupported("2D and solid sections in one model".into())); } let material = model @@ -206,8 +271,17 @@ impl Mesh { .ok_or_else(|| Error::Invalid(format!("section `{name}`: material `{}` is missing", section.material)))?; let elastic = material.elastic.as_ref().ok_or_else(|| Error::Invalid(format!("material `{}` has no elastic constants", section.material)))?; - let d = elasticity(elastic).map_err(|e| Error::Invalid(format!("material `{}`: {e}", section.material)))?; let set = resolved.sets.get(§ion.set).ok_or_else(|| Error::Invalid(format!("set `{}` is not resolved", section.set)))?; + if let Some(area) = rod { + let young = + young_modulus(elastic).map_err(|e| Error::Invalid(format!("material `{}` of section `{name}`: {e}", section.material)))?; + for &label in &set.elements { + let el = model.element(label).ok_or_else(|| Error::Invalid(format!("element {label} is missing")))?; + trusses.push((label, el.kind, el.nodes.clone(), young, area, material.density)); + } + continue; + } + let d = elasticity(elastic).map_err(|e| Error::Invalid(format!("material `{}`: {e}", section.material)))?; for &label in &set.elements { let el = model.element(label).ok_or_else(|| Error::Invalid(format!("element {label} is missing")))?; let (formulation, dofs) = match theory { @@ -222,10 +296,19 @@ impl Mesh { return Err(Error::NotSupported(format!("couplings in 2D models (coupling `{name}`)"))); } elements.sort_by_key(|e| e.0); + trusses.sort_by_key(|e| e.0); + let translations = translation_kinds(planar.is_some()); + // What discrete elements need: at mesh nodes (some only they use), and at reference + // points, which then get their six degrees of freedom like coupled ones. + let needs: Vec<(u32, DofSet)> = crate::discrete::needs(resolved, translations).collect(); let mut labels: Vec = elements.iter().flat_map(|e| e.2.iter().copied()).collect(); + labels.extend(trusses.iter().flat_map(|e| e.2.iter().copied())); + labels.extend(needs.iter().map(|n| n.0).filter(|l| model.node(*l).is_some())); labels.sort_unstable(); labels.dedup(); - let coupled: BTreeSet<&str> = model.couplings().values().map(|c| c.reference_point.as_str()).collect(); + let mut coupled: BTreeSet<&str> = model.couplings().values().map(|c| c.reference_point.as_str()).collect(); + let discrete_labels: BTreeSet = needs.iter().map(|n| n.0).collect(); + coupled.extend(model.reference_points().iter().filter(|(_, p)| discrete_labels.contains(&p.label)).map(|(n, _)| n.as_str())); let points: Vec = coupled .into_iter() .map(|name| { @@ -237,6 +320,12 @@ impl Mesh { for e in &elements { layout.element(&e.2, e.7); } + for e in &trusses { + layout.element(&e.2, modelcraft_fem::ElementDofs::new(translations)); + } + for (label, kinds) in needs.iter().filter(|n| model.node(n.0).is_some()) { + layout.node(*label, *kinds); + } for p in &points { layout.point(p.label, DofSet::STRUCTURAL); } @@ -249,14 +338,46 @@ impl Mesh { .iter() .map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| Error::Invalid(format!("node {l} is missing")))) .collect::>()?; - let mut mesh = Mesh { labels, xyz, elements: Vec::with_capacity(elements.len()), points, layout, theory: planar }; + let mut mesh = Mesh { + labels, + xyz, + elements: Vec::with_capacity(elements.len()), + points, + layout, + theory: planar, + trusses: Vec::with_capacity(trusses.len()), + discrete: crate::discrete::Discrete::default(), + }; for (label, kind, nodes, formulation, d, density, theory, _) in elements { let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); mesh.elements.push(Elem { label, kind, nodes: idx, xyz: exyz, formulation, d, density, theory }); } + for (label, kind, nodes, young, area, density) in trusses { + let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); + let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); + mesh.trusses.push(Bar { label, kind, nodes: idx, xyz: exyz, young, area, density }); + } + mesh.discrete = crate::discrete::Discrete::new(&mesh, resolved)?; Ok(mesh) } + + /// The layout slots (mesh nodes, reference points) every element and discrete element + /// couples, for the matrix pattern. + pub fn blocks(&self) -> impl Iterator> { + let elements = self.elements.iter().map(|e| e.nodes.as_slice()).chain(self.trusses.iter().map(|t| t.nodes.as_slice())); + elements.map(std::borrow::Cow::Borrowed).chain(self.discrete.couplings().map(|dofs| { + let mut slots: Vec = dofs.iter().filter_map(|d| self.layout.locate(*d).map(|(slot, _)| slot)).collect(); + slots.dedup(); + std::borrow::Cow::Owned(slots) + })) + } + + /// The slots of the matrix pattern: the mesh nodes, and the reference points too when a + /// discrete element acts on them (couplings tie points through constraints, not the matrix). + pub fn block_count(&self) -> usize { + if self.discrete.couplings().flatten().any(|d| d >= self.mesh_dofs()) { self.layout.slots() } else { self.labels.len() } + } } /// The prescribed value of every degree of freedom (`None` = free) under a step's boundary @@ -279,8 +400,9 @@ pub fn prescribed( if let Target::ReferencePoint(rp) = bc.target() && mesh.point(rp).is_none() { - warnings - .push(format!("boundary condition `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh; it holds nothing")); + warnings.push(format!( + "boundary condition `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh and no spring or point mass uses; it holds nothing" + )); continue; } for slot in resolved.target_nodes(bc.target()).iter().filter_map(|l| mesh.layout.slot(*l)) { @@ -365,7 +487,9 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ if let Some(rp) = load.reference_point() && mesh.point(rp).is_none() { - return Err(Error::Invalid(format!("load `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh"))); + return Err(Error::Invalid(format!( + "load `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh and no spring or point mass uses" + ))); } // Moments act on the rotations, which only reference points have. let rotation = matches!(load, Load::Moment { .. }); @@ -433,15 +557,36 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ let rho = e.density.ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} has no density", e.label)))?; body(&mut f, mesh, e, acceleration.map(|a| rho * a))?; } + for t in &mesh.trusses { + if chosen.is_some_and(|c| c.binary_search(&t.label).is_err()) { + continue; + } + let rho = t.density.ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} has no density", t.label)))?; + bar_body(&mut f, mesh, t, acceleration.map(|a| rho * a))?; + } + // Point masses take gravity on the whole model. + if set.is_none() { + mesh.discrete.add_gravity(&mut f, *acceleration, mesh); + } } Load::BodyForce { set, vector } => { for label in resolved.sets.get(set).map(|r| r.elements.as_slice()).unwrap_or_default() { + if let Some(t) = mesh.truss(*label) { + bar_body(&mut f, mesh, t, *vector)?; + continue; + } let e = mesh.element(*label).ok_or_else(|| Error::Invalid(format!("load `{name}`: element {label} is not a solid element")))?; body(&mut f, mesh, e, *vector)?; } } Load::Centrifugal { angular_velocity, point, axis, set } => { let chosen: Option<&[u32]> = set.as_ref().map(|s| resolved.sets.get(s).map(|r| r.elements.as_slice()).unwrap_or_default()); + if let Some(t) = mesh.trusses.iter().find(|t| chosen.is_none_or(|c| c.binary_search(&t.label).is_ok())) { + return Err(Error::NotSupported(format!("centrifugal loads on trusses (load `{name}`, element {})", t.label))); + } + if set.is_none() && !mesh.discrete.masses.is_empty() { + return Err(Error::NotSupported(format!("centrifugal loads on point masses (load `{name}`)"))); + } for e in &mesh.elements { if chosen.is_some_and(|c| c.binary_search(&e.label).is_err()) { continue; @@ -461,6 +606,15 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ Ok(f) } +/// A force `b` per volume on a truss. +fn bar_body(f: &mut [f64], mesh: &Mesh, t: &Bar, b: [f64; 3]) -> Result<()> { + let forces = t.truss()?.body_force(b, t.area).map_err(|err| Error::element(t.label, err))?; + for (node, v) in t.nodes.iter().zip(forces) { + add(f, &translations(mesh, *node)?, v); + } + Ok(()) +} + fn body(f: &mut [f64], mesh: &Mesh, e: &Elem, b: [f64; 3]) -> Result<()> { let forces = e.kernel()?.body_force(b).map_err(|err| Error::element(e.label, err))?; for (node, v) in e.nodes.iter().zip(forces) { diff --git a/crates/solver/tests/one_d.rs b/crates/solver/tests/one_d.rs new file mode 100644 index 0000000..ac099b0 --- /dev/null +++ b/crates/solver/tests/one_d.rs @@ -0,0 +1,652 @@ +//! Trusses, point masses, rotary inertia and springs against closed-form solutions. Run with +//! `--nocapture` to see the errors. +#![allow(clippy::needless_range_loop)] +//! +//! - Statically determinate trusses (a plane Pratt truss and a space tripod): the member forces +//! are those of the method of joints, whatever the members' stiffness. +//! - The three-bar truss (statically indeterminate): N = P cos²θ / (1 + 2 cos³θ) in the inclined +//! bars, P / (1 + 2 cos³θ) in the vertical one. +//! - rod3: constant-stress patch with off-centre middle nodes; a bar hanging under its own +//! weight (quadratic displacement, exact at the nodes); curved members are checked in +//! `modelcraft_fem::truss` (a ring's radial expansion gives the hoop strain δ/R exactly). +//! - Spring–mass ω = √(k/m), springs along three skew directions, the two-mass chain, a bar with +//! a tip mass ω = √(EA/(LM)), the torsional pendulum ω = √(k_t/J). + +use std::f64::consts::PI; + +use modelcraft_fem::truss::axial; +use modelcraft_model::*; +use modelcraft_solver::{FrequencyOptions, FrequencyResult, Limits, MassKind, StaticResult, solve_frequency, solve_static}; + +const ZERO_Z: [Option; 6] = [None, None, Some(0.0), None, None, None]; + +/// A model built point by point. +struct Build { + m: Model, + step: &'static str, +} + +impl Build { + fn new(procedure: Procedure) -> Build { + let mut m = Model::default(); + m.add_step(Step::new("Step", procedure), None).unwrap(); + Build { m, step: "Step" } + } + + fn statics() -> Build { + Build::new(Procedure::Static {}) + } + + fn modes(n: u32) -> Build { + Build::new(Procedure::Frequency { modes: n, min_hz: None, max_hz: None }) + } + + fn node(&mut self, label: u32, xyz: [f64; 3]) { + self.m.set_node(label, xyz).unwrap(); + } + + fn nodes(&mut self, name: &str, nodes: &[u32]) -> Target { + self.m.add_set(name, Set { nodes: nodes.to_vec(), ..Set::default() }).unwrap(); + Target::Set(name.into()) + } + + /// Members with one area and material. + fn members(&mut self, name: &str, kind: ElementKind, members: &[(u32, Vec)], area: f64, young: f64, density: Option) { + for (label, nodes) in members { + self.m.set_element(*label, kind, nodes.clone()).unwrap(); + } + self.m.add_set(name, Set { elements: members.iter().map(|m| m.0).collect(), ..Set::default() }).unwrap(); + self.m.add_material(name, Material { elastic: Some(Elastic::Isotropic { young, poisson: 0.3 }), density, ..Material::default() }).unwrap(); + self.m + .add_section( + name, + Section { set: name.into(), material: name.into(), kind: SectionKind::Rod { area }, formulation: Formulation::default() }, + ) + .unwrap(); + } + + fn bc(&mut self, name: &str, target: Target, values: [Option; 6]) { + self.m.add_boundary_condition(name, BoundaryCondition::Displacement { target, values }).unwrap(); + self.m.set_boundary_condition_active(self.step, name, true).unwrap(); + } + + fn fix(&mut self, name: &str, target: Target) { + self.m.add_boundary_condition(name, BoundaryCondition::Fixed { target }).unwrap(); + self.m.set_boundary_condition_active(self.step, name, true).unwrap(); + } + + fn load(&mut self, name: &str, load: Load) { + self.m.add_load(name, load).unwrap(); + self.m.set_load_active(self.step, name, true).unwrap(); + } + + fn force(&mut self, name: &str, node: u32, vector: [f64; 3]) { + let target = self.nodes(name, &[node]); + self.load(name, Load::Force { target, vector, distribution: Distribution::Total }); + } + + fn spring(&mut self, name: &str, first: Target, second: Option, direction: Direction, stiffness: f64) { + self.m.add_spring(name, Spring::Spring { first, second, direction, stiffness }).unwrap(); + } + + fn mass(&mut self, name: &str, target: Target, mass: Option, inertia: Option<[f64; 6]>) { + self.m.add_mass(name, PointMass { target, mass, inertia }).unwrap(); + } + + fn solve(&self) -> StaticResult { + solve_static(&self.m, None, &Limits::default()).unwrap() + } + + fn modes_of(&self, mass: MassKind) -> FrequencyResult { + solve_frequency(&self.m, None, &FrequencyOptions { mass, shift: None }, &Limits::default()).unwrap() + } +} + +fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[0] - b[0], a[1] - b[1], a[2] - b[2]] +} + +fn unit(a: [f64; 3]) -> [f64; 3] { + let l = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt(); + a.map(|c| c / l) +} + +/// The axial force of each member at its first integration point: N = A tᵀ σ t. +fn member_forces(m: &Model, r: &StaticResult, area_of: impl Fn(u32) -> f64) -> Vec<(u32, f64)> { + r.elements + .iter() + .map(|e| { + let el = m.element(e.label).unwrap(); + let t = unit(sub(m.node(el.nodes[1]).unwrap().xyz, m.node(el.nodes[0]).unwrap().xyz)); + (e.label, area_of(e.label) * axial(&e.points[0].stress, t, false)) + }) + .collect() +} + +/// Gaussian elimination with partial pivoting (the method of joints' equations). +fn solve_dense(mut a: Vec>, mut b: Vec) -> Vec { + let n = b.len(); + for c in 0..n { + let p = (c..n).max_by(|i, j| a[*i][c].abs().total_cmp(&a[*j][c].abs())).unwrap(); + a.swap(c, p); + b.swap(c, p); + assert!(a[c][c].abs() > 1e-12, "the truss is not statically determinate"); + for r in c + 1..n { + let f = a[r][c] / a[c][c]; + for k in c..n { + a[r][k] -= f * a[c][k]; + } + b[r] -= f * b[c]; + } + } + let mut x = vec![0.0; n]; + for r in (0..n).rev() { + x[r] = (b[r] - (r + 1..n).map(|k| a[r][k] * x[k]).sum::()) / a[r][r]; + } + x +} + +/// Member forces by the method of joints: at every joint and direction, Σ N t + F + R = 0, with +/// the reactions `supports` (node, direction) as unknowns too. +fn method_of_joints( + xyz: &[(u32, [f64; 3])], + members: &[(u32, [u32; 2])], + supports: &[(u32, usize)], + loads: &[(u32, [f64; 3])], + dims: usize, +) -> Vec<(u32, f64)> { + let row = |node: u32, d: usize| xyz.iter().position(|(l, _)| *l == node).unwrap() * dims + d; + let pos = |node: u32| xyz.iter().find(|(l, _)| *l == node).unwrap().1; + let n = xyz.len() * dims; + assert_eq!(members.len() + supports.len(), n, "determinate"); + let mut a = vec![vec![0.0; n]; n]; + let mut b = vec![0.0; n]; + for (k, (_, [i, j])) in members.iter().enumerate() { + let t = unit(sub(pos(*j), pos(*i))); + for d in 0..dims { + // Tension pulls joint i towards j and j towards i. + a[row(*i, d)][k] += t[d]; + a[row(*j, d)][k] -= t[d]; + } + } + for (k, (node, d)) in supports.iter().enumerate() { + a[row(*node, *d)][members.len() + k] = 1.0; + } + for (node, f) in loads { + for d in 0..dims { + b[row(*node, d)] -= f[d]; + } + } + let x = solve_dense(a, b); + members.iter().zip(x).map(|((l, _), n)| (*l, n)).collect() +} + +fn worst_relative(fe: &[(u32, f64)], exact: &[(u32, f64)]) -> f64 { + let scale = exact.iter().fold(0.0f64, |m, (_, n)| m.max(n.abs())); + exact.iter().map(|(l, n)| (fe.iter().find(|(m, _)| m == l).unwrap().1 - n).abs() / scale).fold(0.0, f64::max) +} + +/// A plane Pratt truss of four panels, pinned at one end and on a roller at the other, loaded at +/// the bottom joints; chords and web members have different areas. The out-of-plane direction is +/// held. Member forces equal the method of joints'. +#[test] +fn pratt_truss_member_forces_by_the_method_of_joints() { + let (a, h, p) = (3.0, 4.0, 10.0); + let mut b = Build::statics(); + let mut xyz = Vec::new(); + for i in 0..5u32 { + xyz.push((1 + i, [a * f64::from(i), 0.0, 0.0])); + } + for i in 1..4u32 { + xyz.push((10 + i, [a * f64::from(i), h, 0.0])); + } + for (l, x) in &xyz { + b.node(*l, *x); + } + let chords: Vec<(u32, [u32; 2])> = vec![(1, [1, 2]), (2, [2, 3]), (3, [3, 4]), (4, [4, 5]), (5, [11, 12]), (6, [12, 13])]; + let web: Vec<(u32, [u32; 2])> = vec![(7, [2, 11]), (8, [3, 12]), (9, [4, 13]), (10, [1, 11]), (11, [13, 5]), (12, [11, 3]), (13, [13, 3])]; + let as_members = |list: &[(u32, [u32; 2])]| list.iter().map(|(l, n)| (*l, n.to_vec())).collect::>(); + b.members("chords", ElementKind::Rod2, &as_members(&chords), 2.0, 2.0e5, None); + b.members("web", ElementKind::Rod2, &as_members(&web), 0.7, 7.0e4, None); + let all: Vec = xyz.iter().map(|x| x.0).collect(); + let plane = b.nodes("all_nodes", &all); + b.bc("plane", plane, ZERO_Z); + let pin = b.nodes("pin", &[1]); + b.bc("pin", pin, [Some(0.0), Some(0.0), None, None, None, None]); + let roller = b.nodes("roller", &[5]); + b.bc("roller", roller, [None, Some(0.0), None, None, None, None]); + let loads = [(2, [0.0, -p, 0.0]), (3, [0.0, -1.5 * p, 0.0]), (4, [0.2 * p, -p, 0.0])]; + for (node, f) in loads { + b.force(&format!("P{node}"), node, f); + } + let r = b.solve(); + let fe = member_forces(&b.m, &r, |l| if l <= 6 { 2.0 } else { 0.7 }); + let members: Vec<(u32, [u32; 2])> = chords.iter().chain(&web).copied().collect(); + let exact = method_of_joints(&xyz, &members, &[(1, 0), (1, 1), (5, 1)], &loads, 2); + let err = worst_relative(&fe, &exact); + eprintln!("Pratt truss: worst member force error {err:.2e} of the largest"); + assert!(err < 1e-10, "{fe:?} vs {exact:?}"); + // Reactions balance the loads. + for d in 0..2 { + assert!((r.applied_total[d] + r.reaction_total[d]).abs() < 1e-9 * p); + } +} + +/// A space tripod: three members from an apex to three fixed points, a skew load at the apex. +#[test] +fn space_tripod_member_forces() { + let (rad, h) = (2.0, 3.0); + let mut b = Build::statics(); + let apex = (4, [0.1, -0.2, h]); + let feet: Vec<(u32, [f64; 3])> = (0..3) + .map(|i| (1 + i as u32, [rad * (PI / 2.0 + 2.0 * PI * i as f64 / 3.0).cos(), rad * (PI / 2.0 + 2.0 * PI * i as f64 / 3.0).sin(), 0.0])) + .collect(); + for (l, x) in feet.iter().chain([&apex]) { + b.node(*l, *x); + } + let members: Vec<(u32, Vec)> = (1..=3).map(|i| (i, vec![i, 4])).collect(); + b.members("legs", ElementKind::Rod2, &members, 1.5, 1.0e5, None); + let feet_set = b.nodes("feet", &[1, 2, 3]); + b.fix("feet", feet_set); + let load = [3.0, -2.0, -10.0]; + b.force("P", 4, load); + let r = b.solve(); + let fe = member_forces(&b.m, &r, |_| 1.5); + // Equilibrium of the apex: Σ N_i t_i + P = 0 with t_i from the apex to foot i (tension + // pulls the apex towards the foot). + let t: Vec<[f64; 3]> = feet.iter().map(|(_, x)| unit(sub(*x, apex.1))).collect(); + let a = (0..3).map(|d| (0..3).map(|i| t[i][d]).collect()).collect(); + let exact: Vec<(u32, f64)> = (1..=3).zip(solve_dense(a, load.map(|c| -c).to_vec())).collect(); + let err = worst_relative(&fe, &exact); + eprintln!("space tripod: worst member force error {err:.2e}"); + assert!(err < 1e-10, "{fe:?} vs {exact:?}"); +} + +/// The statically indeterminate three-bar truss: a vertical bar of length L and two at ±θ +/// meeting at a loaded joint. +#[test] +fn three_bar_truss() { + let (l, theta, p, ea) = (2.0, 35f64.to_radians(), 50.0, 3.0e4); + for kind in [ElementKind::Rod2, ElementKind::Rod3] { + let mut b = Build::statics(); + b.node(1, [-l * theta.tan(), l, 0.0]); + b.node(2, [0.0, l, 0.0]); + b.node(3, [l * theta.tan(), l, 0.0]); + b.node(4, [0.0, 0.0, 0.0]); + let mut members = Vec::new(); + for (i, top) in [1u32, 2, 3].into_iter().enumerate() { + let mut nodes = vec![top, 4]; + if kind == ElementKind::Rod3 { + // Middle nodes a third of the way along: still straight, not centred. + let (x0, x1) = (b.m.node(top).unwrap().xyz, [0.0; 3]); + let mid = 10 + i as u32; + b.node(mid, [0, 1, 2].map(|k| x0[k] + (x1[k] - x0[k]) / 3.0)); + nodes.push(mid); + // A straight member's middle node has no stiffness across it: hold it there. + let across = unit([x0[1], -x0[0], 0.0]); + let set = b.nodes(&format!("mid{mid}"), &[mid]); + b.spring(&format!("mid{mid}"), set, None, Direction::Vector(across), 1.0); + } + members.push((1 + i as u32, nodes)); + } + b.members("bars", kind, &members, 1.0, ea, None); + let top = b.nodes("top", &[1, 2, 3]); + b.fix("top", top); + let all: Vec = b.m.nodes().keys().copied().collect(); + let plane = b.nodes("plane_nodes", &all); + b.bc("plane", plane, ZERO_Z); + b.force("P", 4, [0.0, -p, 0.0]); + let r = b.solve(); + let c3 = theta.cos().powi(3); + let vertical = p / (1.0 + 2.0 * c3); + let inclined = p * theta.cos().powi(2) / (1.0 + 2.0 * c3); + let fe = member_forces(&b.m, &r, |_| 1.0); + let exact = [(1, inclined), (2, vertical), (3, inclined)]; + let err = worst_relative(&fe, &exact); + let delta = r.displacement[r.nodes.binary_search(&4).unwrap()][1]; + let want = -vertical * l / ea; + eprintln!("three-bar truss ({kind:?}): force error {err:.2e}, deflection {delta:.6e} vs {want:.6e}"); + assert!(err < 1e-10); + assert!((delta - want).abs() < 1e-10 * want.abs()); + } +} + +/// A rod3 bar with off-centre middle nodes under end tension carries σ = F/A exactly, at every +/// integration point and extrapolated to every node. +#[test] +fn rod3_constant_stress_patch() { + let (area, young, f) = (2.0, 1000.0, 6.0); + let mut b = Build::statics(); + let dir = unit([1.0, 2.0, 2.0]); + let at = |s: f64| dir.map(|c| c * s); + // Three elements of unequal length, middle nodes at 30%, 50% and 65% of each. + let ends = [0.0, 1.0, 2.5, 3.2]; + let fractions = [0.3, 0.5, 0.65]; + for (i, s) in ends.iter().enumerate() { + b.node(1 + i as u32, at(*s)); + } + let mut members = Vec::new(); + for i in 0..3 { + let s = ends[i] + fractions[i] * (ends[i + 1] - ends[i]); + b.node(10 + i as u32, at(s)); + members.push((1 + i as u32, vec![1 + i as u32, 2 + i as u32, 10 + i as u32])); + } + b.members("bar", ElementKind::Rod3, &members, area, young, None); + let start = b.nodes("start", &[1]); + b.fix("start", start); + // The other nodes slide along the bar only: hold them with springs across it. + let others: Vec = b.m.nodes().keys().copied().filter(|l| *l != 1).collect(); + let rest = b.nodes("rest", &others); + let across = unit([2.0, -1.0, 0.0]); + let third = [dir[1] * across[2] - dir[2] * across[1], dir[2] * across[0] - dir[0] * across[2], dir[0] * across[1] - dir[1] * across[0]]; + b.spring("across1", rest.clone(), None, Direction::Vector(across), 1.0e6); + b.spring("across2", rest, None, Direction::Vector(third), 1.0e6); + b.force("pull", 4, dir.map(|c| c * f)); + let r = b.solve(); + let mut worst: f64 = 0.0; + for e in &r.elements { + for p in &e.points { + worst = worst.max((axial(&p.stress, dir, false) - f / area).abs()); + } + for s in &e.nodal { + worst = worst.max((axial(s, dir, false) - f / area).abs()); + } + } + let tip = r.displacement[r.nodes.binary_search(&4).unwrap()]; + let stretch = tip.iter().zip(dir).map(|(u, d)| u * d).sum::(); + eprintln!("rod3 patch: stress error {worst:.2e}; tip stretch {stretch:.10} vs FL/EA {:.10}", f * 3.2 / (young * area)); + assert!(worst < 1e-10 * f / area); + assert!((stretch - f * 3.2 / (young * area)).abs() < 1e-12); +} + +/// A bar hanging under its own weight: u(z) = ρ g (L z − z²/2) / E from the support. Two-node +/// bars are exact at the nodes in one dimension; three-node bars reproduce the quadratic exactly. +#[test] +fn hanging_bar_under_gravity() { + let (l, rho, g, young, area) = (6.0, 2.0, 9.81, 1.0e4, 0.5); + for (kind, cells) in [(ElementKind::Rod2, 6usize), (ElementKind::Rod3, 3)] { + let mut b = Build::statics(); + let h = l / cells as f64; + for i in 0..=cells { + b.node(1 + i as u32, [0.0, 0.0, -(i as f64) * h]); + } + let mut members = Vec::new(); + for i in 0..cells { + let mut nodes = vec![1 + i as u32, 2 + i as u32]; + if kind == ElementKind::Rod3 { + b.node(100 + i as u32, [0.0, 0.0, -(i as f64 + 0.5) * h]); + nodes.push(100 + i as u32); + } + members.push((1 + i as u32, nodes)); + } + b.members("bar", kind, &members, area, young, Some(rho)); + let top = b.nodes("top", &[1]); + b.fix("top", top); + let all: Vec = b.m.nodes().keys().copied().collect(); + let line = b.nodes("line", &all); + b.bc("line", line, [Some(0.0), Some(0.0), None, None, None, None]); + b.load("weight", Load::Gravity { acceleration: [0.0, 0.0, -g], set: None }); + let r = b.solve(); + let mut worst: f64 = 0.0; + for (label, u) in r.nodes.iter().zip(&r.displacement) { + let z = -b.m.node(*label).unwrap().xyz[2]; + let exact = -rho * g * (l * z - z * z / 2.0) / young; + worst = worst.max((u[2] - exact).abs()); + } + let tip = rho * g * l * l / (2.0 * young); + eprintln!("hanging bar ({kind:?}): nodal displacement error {:.2e} of the tip's", worst / tip); + assert!(worst < 1e-12 * tip.max(1.0), "{worst}"); + assert!((r.reaction_total[2] - rho * g * area * l).abs() < 1e-10); + } +} + +/// One node, a spring to the ground in X and a mass: ω = √(k/m) (consistent and lumped mass +/// alike); a static force gives u = F/k. +#[test] +fn spring_mass_natural_frequency_and_static_stretch() { + let (k, m) = (400.0, 2.5); + for mass in [MassKind::Consistent, MassKind::Lumped] { + let mut b = Build::modes(1); + b.node(1, [1.0, 2.0, 3.0]); + let n = b.nodes("n", &[1]); + b.spring("k", n.clone(), None, Direction::Dof(1), k); + b.mass("m", n.clone(), Some(Mass::Isotropic(m)), None); + b.bc("yz", n, [None, Some(0.0), Some(0.0), None, None, None]); + let r = b.modes_of(mass); + let f = (k / m).sqrt() / (2.0 * PI); + let err = (r.modes[0].frequency - f).abs() / f; + eprintln!("spring–mass ({mass:?}): f {:.12} vs √(k/m)/2π {f:.12}, error {err:.1e}", r.modes[0].frequency); + assert!(err < 1e-12); + assert!((r.total_mass - m).abs() < 1e-12); + assert!((r.modes[0].effective_mass[0] - m).abs() < 1e-9); + } + let mut b = Build::statics(); + b.node(1, [0.0; 3]); + let n = b.nodes("n", &[1]); + b.spring("k", n.clone(), None, Direction::Dof(1), k); + b.bc("yz", n, [None, Some(0.0), Some(0.0), None, None, None]); + b.force("F", 1, [8.0, 0.0, 0.0]); + let r = b.solve(); + assert!((r.displacement[0][0] - 8.0 / k).abs() < 1e-15); + // The spring carries the load to the ground: no support reaction, and no imbalance. + assert_eq!(r.reaction_total[0], 0.0); + assert!(r.warnings.iter().all(|w| !w.contains("differ")), "{:?}", r.warnings); +} + +/// Springs to the ground along three orthogonal skew vectors: frequencies √(kᵢ/m) with the mode +/// shapes along the vectors. +#[test] +fn springs_in_skew_directions() { + let m = 3.0; + let a = unit([1.0, 1.0, 0.5]); + let b2 = unit([1.0, -1.0, 0.0]); + let c = [a[1] * b2[2] - a[2] * b2[1], a[2] * b2[0] - a[0] * b2[2], a[0] * b2[1] - a[1] * b2[0]]; + let ks = [30.0, 120.0, 270.0]; + let mut b = Build::modes(3); + b.node(7, [0.0; 3]); + let n = b.nodes("n", &[7]); + for (i, (v, k)) in [a, b2, c].iter().zip(ks).enumerate() { + b.spring(&format!("k{i}"), n.clone(), None, Direction::Vector(v.map(|x| 2.0 * x)), k); + } + b.mass("m", n, Some(Mass::Isotropic(m)), None); + let r = b.modes_of(MassKind::Consistent); + for ((mode, k), v) in r.modes.iter().zip(ks).zip([a, b2, c]) { + let f = (k / m).sqrt() / (2.0 * PI); + let s = mode.shape[0]; + let along = (s[0] * v[0] + s[1] * v[1] + s[2] * v[2]).abs() / (s[0] * s[0] + s[1] * s[1] + s[2] * s[2]).sqrt(); + eprintln!("skew springs: f {:.10} vs {f:.10}; shape along the spring {along:.12}", mode.frequency); + assert!((mode.frequency - f).abs() < 1e-10 * f); + assert!((along - 1.0).abs() < 1e-9); + } +} + +/// Ground — k₁ — m₁ — k₂ — m₂, the two-point spring acting along the line of its nodes: +/// ω² are the roots of m₁m₂ω⁴ − (k₁m₂ + k₂m₁ + k₂m₂)ω² + k₁k₂ = 0. +#[test] +fn two_mass_chain() { + let (k1, k2, m1, m2) = (100.0, 60.0, 2.0, 0.5); + let mut b = Build::modes(2); + let dir = unit([3.0, 0.0, 4.0]); + b.node(1, [0.0; 3]); + b.node(2, dir.map(|c| 2.0 * c)); + let n1 = b.nodes("n1", &[1]); + let n2 = b.nodes("n2", &[2]); + b.spring("k1", n1.clone(), None, Direction::Vector(dir), k1); + b.spring("k2", n1.clone(), Some(n2.clone()), Direction::Line, k2); + b.mass("m1", n1, Some(Mass::Isotropic(m1)), None); + b.mass("m2", n2, Some(Mass::Directional([m2, 7.0, m2])), None); + // Only the motion along the line: hold the others with stiff springs across it. + let both = b.nodes("both", &[1, 2]); + b.spring("hold_y", both.clone(), None, Direction::Dof(2), 1.0e9); + b.spring("hold_n", both, None, Direction::Vector([dir[2], 0.0, -dir[0]]), 1.0e9); + let r = b.modes_of(MassKind::Lumped); + let (pa, pb, pc) = (m1 * m2, -(k1 * m2 + k2 * m1 + k2 * m2), k1 * k2); + let disc = (pb * pb - 4.0 * pa * pc).sqrt(); + let exact = [(-pb - disc) / (2.0 * pa), (-pb + disc) / (2.0 * pa)]; + for (mode, w2) in r.modes.iter().zip(exact) { + let err = (mode.eigenvalue - w2).abs() / w2; + eprintln!("two-mass chain: ω² {:.10} vs {w2:.10}, error {err:.1e}", mode.eigenvalue); + assert!(err < 1e-8, "the holding springs perturb at 1e-9"); + } +} + +/// A massless bar with a tip mass M: ω = √(EA/(L M)); the bar's own mass ρAL adds the +/// Rayleigh-quotient correction ω² ≈ EA/(L(M + ρAL/3)) for a light bar. +#[test] +fn bar_with_a_tip_mass() { + let (l, young, area, tip) = (2.0, 2.0e5, 1e-2, 50.0); + for kind in [ElementKind::Rod2, ElementKind::Rod3] { + let mut b = Build::modes(1); + b.node(1, [0.0; 3]); + b.node(2, [l, 0.0, 0.0]); + let mut nodes = vec![1, 2]; + if kind == ElementKind::Rod3 { + b.node(3, [l / 2.0, 0.0, 0.0]); + nodes.push(3); + } + b.members("bar", kind, &[(1, nodes)], area, young, Some(1e-9)); + let root = b.nodes("root", &[1]); + b.fix("root", root); + let all: Vec = b.m.nodes().keys().copied().collect(); + let line = b.nodes("line", &all); + b.bc("line", line, [None, Some(0.0), Some(0.0), None, None, None]); + let end = b.nodes("end", &[2]); + b.mass("tip", end, Some(Mass::Isotropic(tip)), None); + let r = b.modes_of(MassKind::Consistent); + let f = (young * area / (l * tip)).sqrt() / (2.0 * PI); + let err = (r.modes[0].frequency - f).abs() / f; + eprintln!("bar with tip mass ({kind:?}): f {:.10} vs {f:.10}, error {err:.1e}", r.modes[0].frequency); + assert!(err < 1e-9); + } +} + +/// A torsional pendulum: a reference point with polar inertia J on a torsional spring k_t to +/// the ground; ω = √(k_t/J). The other directions are held. +#[test] +fn torsional_pendulum() { + let (kt, j) = (5.0e3, 0.8); + let mut b = Build::modes(1); + b.m.add_reference_point("Disc", ReferencePoint { xyz: [1.0, 1.0, 0.0], label: 0 }).unwrap(); + let disc = Target::ReferencePoint("Disc".into()); + b.spring("shaft", disc.clone(), None, Direction::Dof(6), kt); + // Products of inertia about the other axes do not couple into the held rotations. + b.mass("disc", disc.clone(), Some(Mass::Isotropic(4.0)), Some([0.5, 0.5, j, 0.0, 0.1, 0.1])); + b.bc("held", disc, [Some(0.0), Some(0.0), Some(0.0), Some(0.0), Some(0.0), None]); + let r = b.modes_of(MassKind::Lumped); + let f = (kt / j).sqrt() / (2.0 * PI); + let err = (r.modes[0].frequency - f).abs() / f; + eprintln!("torsional pendulum: f {:.12} vs √(k_t/J)/2π {f:.12}, error {err:.1e}", r.modes[0].frequency); + assert!(err < 1e-12); + assert_eq!(r.reference_points.len(), 1); + assert!(r.modes[0].point_shape[0][5].abs() > 0.0); +} + +/// Gravity acts on point masses (whole-model gravity); dashpots are ignored with a warning in +/// a static step; the runs are bit-for-bit repeatable. +#[test] +fn point_mass_gravity_dashpots_and_determinism() { + let mut b = Build::statics(); + b.node(1, [0.0; 3]); + b.node(2, [0.0, 0.0, -1.0]); + b.members("rod", ElementKind::Rod2, &[(1, vec![1, 2])], 1.0, 1.0e3, Some(0.0 + 1e-6)); + let top = b.nodes("top", &[1]); + b.fix("top", top); + let bottom = b.nodes("bottom", &[2]); + b.bc("line", bottom.clone(), [Some(0.0), Some(0.0), None, None, None, None]); + b.mass("weight", bottom.clone(), Some(Mass::Isotropic(3.0)), None); + b.m.add_spring("damper", Spring::Dashpot { first: bottom, second: None, direction: Direction::Dof(3), coefficient: 2.0 }).unwrap(); + b.load("g", Load::Gravity { acceleration: [0.0, 0.0, -10.0], set: None }); + let r = b.solve(); + let total = 3.0 * 10.0 + 1e-6 * 10.0; + assert!((r.reaction_total[2] - total).abs() < 1e-10, "{:?}", r.reaction_total); + assert!((r.displacement[r.nodes.binary_search(&2).unwrap()][2] + total * 1.0 / 1.0e3).abs() < 1e-6); + assert!(r.warnings.iter().any(|w| w.contains("dashpots `damper`")), "{:?}", r.warnings); + assert_eq!(r, b.solve()); +} + +/// What the solver cannot do is an error that says so. +#[test] +fn errors_name_the_problem() { + // A point mass with nothing holding it. + let mut b = Build::modes(1); + b.node(1, [0.0; 3]); + let n = b.nodes("n", &[1]); + b.mass("m", n, Some(Mass::Isotropic(1.0)), None); + // A free particle has no stiffness at all to shift against. + assert!(solve_frequency(&b.m, None, &FrequencyOptions::default(), &Limits::default()).is_err()); + // A truss of an orthotropic material. + let mut b = Build::statics(); + b.node(1, [0.0; 3]); + b.node(2, [1.0, 0.0, 0.0]); + b.m.set_element(1, ElementKind::Rod2, vec![1, 2]).unwrap(); + b.m.add_set("rod", Set { elements: vec![1], ..Set::default() }).unwrap(); + b.m.add_material( + "Ortho", + Material { + elastic: Some(Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [1.0; 3] }), + density: None, + ..Material::default() + }, + ) + .unwrap(); + b.m.add_section( + "rod", + Section { set: "rod".into(), material: "Ortho".into(), kind: SectionKind::Rod { area: 1.0 }, formulation: Formulation::default() }, + ) + .unwrap(); + let e = solve_static(&b.m, None, &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("orthotropic"), "{e}"); +} + +/// A 2D (plane stress) plate stiffened by trusses along its loaded edges: the truss takes the +/// model's two in-plane unknowns per node, and the two act as parallel springs, +/// u = F L / (E t w + E_b A) (exact: uniform tension, ν = 0). +#[test] +fn plane_truss_stiffens_a_plane_stress_plate() { + let (l, w, t, young, area, young_bar, force) = (4.0, 1.0, 0.2, 1000.0, 0.05, 2.0e4, 30.0); + for (kind, mid) in [(ElementKind::Rod2, false), (ElementKind::Rod3, true)] { + let mut b = Build::statics(); + for (label, xyz) in [(1, [0.0, 0.0, 0.0]), (2, [l, 0.0, 0.0]), (3, [l, w, 0.0]), (4, [0.0, w, 0.0])] { + b.node(label, xyz); + } + b.m.set_element(1, ElementKind::Quad4, vec![1, 2, 3, 4]).unwrap(); + b.m.add_set("plate", Set { elements: vec![1], ..Set::default() }).unwrap(); + b.m.add_material("Plate", Material { elastic: Some(Elastic::Isotropic { young, poisson: 0.0 }), density: None, ..Material::default() }) + .unwrap(); + b.m.add_section( + "plate", + Section { + set: "plate".into(), + material: "Plate".into(), + kind: SectionKind::PlaneStress { thickness: t }, + formulation: Formulation::default(), + }, + ) + .unwrap(); + // Two bars of half the area, along the bottom and top edges, keep the tension uniform. + let mut members = vec![(2, vec![1, 2]), (3, vec![4, 3])]; + if mid { + b.node(5, [l / 2.0, 0.0, 0.0]); + b.node(6, [l / 2.0, w, 0.0]); + members[0].1.push(5); + members[1].1.push(6); + } + b.members("bar", kind, &members, area / 2.0, young_bar, None); + let left = b.nodes("left", &[1, 4]); + b.bc("left", left, [Some(0.0), None, None, None, None, None]); + let corner = b.nodes("corner", &[1]); + b.bc("corner", corner, [None, Some(0.0), None, None, None, None]); + if mid { + let middle = b.nodes("middle", &[5, 6]); + b.bc("middle", middle, [None, Some(0.0), None, None, None, None]); + } + let right = b.nodes("right", &[2, 3]); + b.load("pull", Load::Force { target: right, vector: [force, 0.0, 0.0], distribution: Distribution::PerNode }); + let r = b.solve(); + let u = r.displacement[r.nodes.binary_search(&2).unwrap()][0]; + let exact = 2.0 * force * l / (young * t * w + young_bar * area); + eprintln!("plane truss ({kind:?}): u {u:.12} vs {exact:.12}; {} unknowns", r.stats.dofs); + assert!((u - exact).abs() < 1e-12 * exact); + assert_eq!(r.stats.dofs, 2 * b.m.nodes().len(), "two unknowns per node in a 2D model"); + } +} diff --git a/crates/ui-egui/src/theme.rs b/crates/ui-egui/src/theme.rs index d56f8f0..55dda3e 100644 --- a/crates/ui-egui/src/theme.rs +++ b/crates/ui-egui/src/theme.rs @@ -153,7 +153,7 @@ impl Tokens { ElementKind::Tet4 | ElementKind::Tet10 => self.tet, ElementKind::Wedge6 => self.wedge, ElementKind::Tri3 | ElementKind::Tri6 | ElementKind::Quad4 | ElementKind::Quad8 => self.shell, - ElementKind::Rod2 | ElementKind::Beam2 => self.line_element, + ElementKind::Rod2 | ElementKind::Rod3 | ElementKind::Beam2 => self.line_element, } } } diff --git a/docs/architecture.md b/docs/architecture.md index 80c0278..ece476f 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -138,8 +138,13 @@ module docs; ADR 0001 records why the linear algebra is our own. quad4, quad8), chosen by the section kind, and the model then has two degrees of freedom per node: the 2D elements declare `ElementDofs::new(DofSet::PLANE)` (Ux, Uy) to the `DofLayout`. A 2D model lies in the XY plane (axisymmetric: x the radius, Y the axis; nodes on the axis are held - radially), is loaded on element edges, and cannot mix with solids or couplings yet. Every reference point a coupling - ties to the mesh gets six degrees of freedom (translations, rotations) after the mesh's. Shells + radially), is loaded on element edges, and cannot mix with solids or couplings yet. Rod + sections' rod2 and rod3 elements are trusses (`modelcraft_fem::truss`). Point masses, rotary + inertia and springs (`Resolved::masses`, `springs`, expanded per node by the model) assemble on + the unknowns their nodes are given (`discrete::needs`, `discrete::dof_of`); dashpots are + ignored with a warning (no velocities in static and frequency steps). Every reference point a + coupling or a discrete element uses gets six degrees of freedom (translations, rotations) after + the mesh's. Shells and beams return "not supported yet". The step's boundary conditions prescribe translations (rotations do not exist on solid nodes) and a coupled point's six directions; its loads become a load vector: forces and tractions from `Resolved::nodal_loads`, forces and moments on coupled @@ -181,7 +186,10 @@ and rigid-body modes. (kinematic and distributing), clamping through a coupled point equals clamping (static and modes), tension through a distributing coupling against the equivalent pressure away from the end, pure bending by a moment through a distributing coupling, torsion of a tube by a moment on -a kinematically coupled point (twist and shear stress). +a kinematically coupled point (twist and shear stress). `crates/solver/tests/one_d.rs`: member +forces of determinate trusses (a plane Pratt truss, a space tripod) against the method of joints, +the three-bar truss, a rod3 constant-stress patch and a bar hanging under its weight, spring–mass, +skew springs, the two-mass chain, a bar with a tip mass and the torsional pendulum. ## Explicit solver @@ -192,9 +200,11 @@ number their degrees of freedom with `modelcraft_fem::DofLayout` (docs/fem.md). 1. **System.** `Model::resolve()`; hex8 sections with reduced integration become one-point uniform-strain hexahedra with Flanagan–Belytschko hourglass control, tet4 sections Bonet–Burton - average nodal pressure tetrahedra (`modelcraft_fem::uniform` holds their kinematics). Other - element kinds and formulations, couplings, reference points and moments return "not supported - yet". Masses are row-sum lumped. A zero-valued boundary condition holds a degree of freedom; a + average nodal pressure tetrahedra (`modelcraft_fem::uniform` holds their kinematics), rod2 and + rod3 trusses (`discrete`). Point masses add to the lumped mass; springs and dashpots (on the + previous half increment's velocity) add internal and viscous forces and their own stable + 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. 2. **Time step.** Each element's highest frequency is bounded in closed form (exact for a cube); @@ -227,7 +237,9 @@ only): stress waves in hex8 and tet4 bars (front at c·t, σ = ρ c v, work inpu rigid wall (contact time 2L/c, rebound speed, plateau force ρ c v₀ A), the energy balance, bit-for-bit repeatability; Dienes' simple shear for both rates, a stressed element turned rigidly 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; in the engine (`solve.rs`), a +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 cantilever's free-vibration period against `solve.frequency`. ## Lattice Boltzmann solver diff --git a/docs/elements.md b/docs/elements.md index 0ede686..9a59649 100644 --- a/docs/elements.md +++ b/docs/elements.md @@ -19,10 +19,10 @@ hex8/20/27 with first and second derivatives). | Element | Use | Status | |---|---|---| -| Point mass | Lumped mass at a node or reference point (isotropic or per direction) | in progress (1D elements I) | -| Rotary inertia | Inertia tensor at a 6-DOF node | in progress (1D elements I) | -| Spring to ground / two-point spring | Linear spring in a fixed direction or along the line of nodes | in progress (1D elements I) | -| Dashpot to ground / two-point | Linear viscous damper (dynamic steps; ignored in static) | in progress (1D elements I) | +| Point mass | Lumped mass at a node or reference point (isotropic or per direction) | done (implicit, explicit; Abaqus MASS, Nastran CONM1, LS-DYNA ELEMENT_MASS; per direction: Nastran only) | +| Rotary inertia | Inertia tensor at a 6-DOF node | done (implicit; Abaqus ROTARYI, Nastran CONM1, LS-DYNA ELEMENT_INERTIA; explicit when it has rotations) | +| Spring to ground / two-point spring | Linear spring in a fixed direction or along the line of nodes | done (implicit, explicit translational; Abaqus SPRING1/2/A, Nastran CELAS2, LS-DYNA ELEMENT_DISCRETE; skew fixed directions: LS-DYNA only) | +| Dashpot to ground / two-point | Linear viscous damper (dynamic steps; ignored in static) | done (explicit; ignored with a warning in static and frequency steps; Abaqus DASHPOT1/2/A, Nastran CDAMP2, LS-DYNA ELEMENT_DISCRETE) | | Nonstructural mass | Mass spread over a region without stiffness | planned | | Connectors | Kinematic joints (hinge, slot, universal…) with elasticity, damping, stops | planned (M7) | @@ -30,8 +30,8 @@ hex8/20/27 with first and second derivatives). | Element | Use | Status | |---|---|---| -| Truss, 2-node (2D and 3D) | Axial members, cables (compression allowed) | in progress (1D elements I) | -| Truss, 3-node (2D and 3D) | Curved axial members | in progress (1D elements I) | +| Truss, 2-node (2D and 3D) | Axial members, cables (compression allowed) | done (rod2: implicit, in plane models with their two unknowns per node; explicit, corotational under nonlinear geometry; Abaqus T3D2 / T2D2, Nastran CROD, LS-DYNA truss beam (3D models)) | +| Truss, 3-node (2D and 3D) | Curved axial members | done (rod3: implicit, plane models too; explicit with small displacements; Abaqus T3D3 / T2D3; no Nastran or LS-DYNA counterpart) | | Beam, 2-node Timoshenko, 3D | Frames, stiffeners; shear deformable | in progress (shells and beams) | | Beam, 2-node Euler–Bernoulli (cubic), 2D and 3D | Slender frames; exact for end-loaded prismatic members | planned (1D elements II) | | Beam, 3-node Timoshenko (quadratic), 2D and 3D | Curved and stocky beams | planned (1D elements II) | diff --git a/docs/features.md b/docs/features.md index 85421da..ac5ce75 100644 --- a/docs/features.md +++ b/docs/features.md @@ -5,7 +5,7 @@ Generated from the capability catalogue (`crates/engine/src/catalogue.rs`) by Every entry is a toolbar button. **Ready** entries run their command; **coming soon** entries are shown greyed with the roadmap milestone that brings them (ROADMAP.md). -48 ready, 53 coming soon. +50 ready, 51 coming soon. ## Geometry @@ -106,8 +106,8 @@ shown greyed with the roadmap milestone that brings them (ROADMAP.md). | Contact | Find Contacts — Detect contact pairs automatically | coming soon (M7) | | | Contact | Rigid Wall — Rigid walls for explicit dynamics | coming soon (M10) | | | Connectors | Bolt Load — Bolt pretension | coming soon (M7) | | -| Connectors | Spring — Springs and dashpots | coming soon (M7) | | -| Connectors | Point Mass — Point masses and inertia | coming soon (M7) | | +| Connectors | Spring — Spring or dashpot to the ground or between two points, along a direction or their line | ready | `model.spring.create`, `model.spring.edit`, `model.spring.delete` | +| Connectors | Point Mass — Point mass and rotary inertia at nodes or a reference point | ready | `model.mass.create`, `model.mass.edit`, `model.mass.delete` | | Connectors | Connector — Connector elements between points | coming soon (M7) | | | Coupled | FSI Interface — The fluid–structure interface | coming soon (M13) | | diff --git a/docs/fem.md b/docs/fem.md index 38d0a86..3e2fb41 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -39,6 +39,12 @@ Lookups: `node(label)` / `slot(label)` (a slot is a position: mesh nodes `0..nod `point_slot(k)`), `dof(slot, kind)`, `range(slot)`, `kinds(slot)`, `locate(dof) -> (slot, kind)`, `element_dofs(slots, kinds)`, `describe(dof)`, `stride()` (uniform layouts). +Trusses declare `modelcraft_fem::truss::DOFS` (translations, on `Basis::Line2`/`Line3`). +Point masses, springs and dashpots are not mesh elements: the implicit solver asks the layout +for what each needs at its nodes (`modelcraft_solver::discrete::needs`: translations for a mass +or a translational spring, rotations for an inertia or a torsional spring; a reference point +they use is a point with six) and finds its unknowns with `discrete::dof_of`. + Solvers: `modelcraft_solver::system::Mesh::layout` numbers the implicit system (mesh nodes, then six per coupled reference point); assembly (`sparse::add_block(k, layout, slots, kinds, data)`, `add_scalar`, `add_diagonal`) works for any mix of kinds per node, so a solid sharing diff --git a/docs/parity.md b/docs/parity.md index 4352e52..54cd3a6 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -56,7 +56,8 @@ products are run, no documentation text is copied). Status: **done**, **started* | Model | Steps: nonlinear static, explicit dynamics, heat transfer; amplitudes | M7–M8 | | planned | | Interaction | Tie between non-matching meshes | M7 | | planned | | Interaction | Contact (general and pairs), friction, automatic pair detection | M7 | | planned | -| Interaction | Bolt pretension, springs, point masses, connectors | M7 | | planned | +| 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 | | Decks | Nastran bulk data writer (`deck.write`, SOL 101 or 103 per deck; reports what Nastran cannot express) | M3 | | done | | Decks | LS-DYNA keyword writer (implicit; explicit in M8) | M3 | | done (implicit static and eigenvalue; beam sections planned) | diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index 68a3537..f87fec7 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -997,6 +997,12 @@ the module docs): distributing (interpolation) constraint as a least-squares rigid fit with loads spread by virtual work, G. H. K. Heirman, W. Desmet, Multibody Syst. Dyn. 24(2):219–234, 2010; the reduced matrix product T_fᵀ K T_f row by row, F. G. Gustavson, ACM TOMS 4(3):250–269, 1978. +- Trusses and discrete elements: the isoparametric bar along its arc, R. D. Cook, D. S. Malkus, + M. E. Plesha, R. J. Witt, *Concepts and Applications of Finite Element Analysis*, 4th ed., + Wiley 2002, §3.6 and §6.2, and Hughes 2000, §1.10; springs k [n nᵀ, −n nᵀ; −n nᵀ, n nᵀ], + Cook et al. 2002, §2.2 and §2.10; HRZ lumping as above. Verified against the method of joints + (statically determinate trusses), the three-bar truss, the hanging bar, √(k/m), the two-mass + chain and the torsional pendulum √(k_t/J). - Verification: R. H. MacNeal, R. L. Harder, *A proposed standard set of problems to test finite element accuracy*, Finite Elem. Anal. Des. 1(1):3–20, 1985 (patch tests); Lamé, Kirsch, Timoshenko and Euler–Bernoulli closed forms; Saint-Venant's principle (tension through a @@ -1024,6 +1030,10 @@ again in the module docs): 14(5):437–449, 1998. - Bulk viscosity and its time-step correction: J. von Neumann, R. D. Richtmyer, J. Appl. Phys. 21(3):232–237, 1950; M. L. Wilkins, J. Comput. Phys. 36(3):281–303, 1980. +- Springs and dashpots: lagged-velocity damping and the damped stable increment + Δt ≤ (2/ω)(√(1 + ξ²) − ξ), Belytschko et al. 2014, §6.2–6.3; Gershgorin's bound on + M^{-1/2} K M^{-1/2} per degree of freedom. Verified against the damped oscillator's decay e^{−ξωt} + and period 2π/ω_d. - Rigid walls as kinematic velocity constraints: T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier, W. Kanoknukulchai, CMAME 8(3):249–276, 1976. - Verification: d'Alembert's solution of the 1D wave equation (front at c t, σ = ρ c v), the