OneCAD-Simulation/docs/fem.md

FEM foundation: degrees of freedom and shape functions

What every element in crates/fem and every solver builds on. Code: crates/fem/src/dof.rs, crates/fem/src/basis.rs, crates/fem/src/jet.rs.

Degrees of freedom

Every unknown is a pair (node, DofKind):

DofKind	Meaning	Model / deck number
Ux, Uy, Uz	translations along X, Y, Z	1, 2, 3
Rx, Ry, Rz	rotations about X, Y, Z	4, 5, 6
P	pressure (mixed u–p elements)	8
T	temperature	11

The enum order (Ux, Uy, Uz, Rx, Ry, Rz, T, P) is the numbering order within a node. DofKind::from_number / number convert to the 1-based numbers boundary conditions use; from_direction / direction to the structural directions 0–5.

Elements declare their unknowns with ElementDofs { node: DofSet, internal: Internal }: the kinds every node of the element carries (DofSet::TRANSLATIONS for solids, PLANE for 2D elements, STRUCTURAL for shells and beams, TEMPERATURE for heat transfer, translations plus PRESSURE for continuous-pressure mixed elements) and the internal unknowns the element condenses out itself (incompatible_modes, discontinuous pressure), which never reach the global system. An element's matrices are ordered node by node, its node kinds in DofKind order within a node (a shell's 24 × 24 stiffness: node 1 Ux Uy Uz Rx Ry Rz, node 2 …). Solid elements: Formulation::dofs() / Solid::dofs().

Numbering (DofLayout): collect what every node needs with DofLayout::builder() (element(labels, dofs), node(label, kinds) for constraints and boundary conditions that need a kind no element gives, point(label, kinds) for reference points), then build(). A node carries the union of what was asked for. Unknowns are numbered mesh nodes in label order, each node's kinds in DofKind order, then the points in the order added; a slot's unknowns are contiguous, so a node pair is one dense block in the assembled matrix. The same model always gives the same numbering.

Lookups: node(label) / slot(label) (a slot is a position: mesh nodes 0..nodes(), points 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 fieldcraft_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 (fieldcraft_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: fieldcraft_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 nodes with a shell adds its translation block into the shell nodes' six-unknown blocks. Boundary conditions and loads address (node, kind): a boundary condition holds the kinds a node carries, and a moment needs rotations. The explicit solver has solid nodes only and checks that its layout is the uniform node-major one (NODE_DOFS = 3).

Shape functions

Basis is every family's interpolation:

Family	Reference element	Nodes
Line2, Line3	ξ ∈ [−1, 1]	ends −1, 1; line3 adds the middle
Tri3, Tri6	ξ, η ≥ 0, ξ + η ≤ 1	corners (0,0), (1,0), (0,1); mid-sides of 1-2, 2-3, 3-1
Quad4, Quad8, Quad9	[−1, 1]²	corners counter-clockwise from (−1,−1); mid-sides of 1-2, 2-3, 3-4, 4-1; centre
Tet4, Tet10	unit tetrahedron	corners (0,0,0), (1,0,0), (0,1,0), (0,0,1); mid-sides of 1-2, 2-3, 3-1, 1-4, 2-4, 3-4
Wedge6, Wedge15	triangle × ζ ∈ [−1, 1]	bottom 1-3, top 4-6; mid-sides of the bottom edges, top edges, vertical edges
Pyramid5	base [−1, 1]² at ζ = 0, apex (0, 0, 1)	base corners as quad4, then the apex
Hex8, Hex20, Hex27	[−1, 1]³	corners 1-4 at ζ = −1, 5-8 above; hex20 mid-sides (bottom, top, vertical edges); hex27 face centres in ElementKind::Hex8.faces() order, then the centre

Node orders are the model's (ElementKind::edges() / faces() are the single source): corners, mid-side nodes in edge order, face nodes, body nodes. Basis::of(ElementKind) maps a model element; nodes() gives natural node coordinates; dim(), count(), corners(), degree() (the complete polynomial degree reproduced: 1 for linear families including the rational pyramid, 2 for the quadratic ones).

Evaluation at a natural point xi: [f64; 3] (unused coordinates are zero):

• eval(xi) -> Shape: values n and first derivatives dn[i] = ∂Nᵢ/∂(ξ, η, ζ). tet4, tet10, wedge6, hex8, hex20 use the hand-written functions of shape.rs (the solid kernels' path).
• eval2(xi) -> Shape2: adds d2n[i][a][b] = ∂²Nᵢ/∂ξₐ∂ξ_b, from second-order forward-mode automatic differentiation of the same formulas (jet::Jet): exact to round-off.
• gradients(&shape, xyz) -> Gradients: physical dx[i] = ∇ₓNᵢ = J⁻¹ ∇_ξNᵢ and det J, with J_ab = ∂x_b/∂ξₐ; an inverted or collapsed element is Error::Degenerate.
• hessians(&shape2, xyz) -> Hessians: physical d2x[i] = J⁻¹ (H_ξ − Σ_b (∂N/∂x_b) X_b″) J⁻ᵀ with X_b″ = Σᵢ x_{i,b} H_ξ(Nᵢ), the map's second derivative (T. J. R. Hughes, The Finite Element Method, 2000, §3.9, §4.2), and laplacian(i). Without the X″ term a field linear in x would have a nonzero Hessian on a curved element.

A family of dimension d maps onto the first d physical coordinates (plane elements in x–y, lines along x). Elements embedded in a higher dimension (shells, beams, trusses in 3D) pass node coordinates in their own local frame.

Verification (basis::tests, dof::tests, jet::tests)

For every family: Kronecker delta at the nodes; partition of unity (Σ N = 1, Σ ∇N = 0, Σ H = 0); exact reproduction of every monomial up to degree() with its gradient and Hessian in natural coordinates; first and second derivatives against central differences; on an affine map, physical gradients, Hessians and Laplacians exact for a complete polynomial of the family's degree; on a curved map, linear fields exact with zero Hessian, and the physical Hessian of a quadratic interpolant consistent with finite differences of its physical gradient. The solver and explicit migration to DofLayout was checked bit for bit: every result of the solver, explicit, engine and io test suites hashes identically before and after.

Heat transfer and thermal strain

A heat conduction element carries one temperature per node (heat::HEAT_DOFS, kind T); heat::Conduction gives the conductivity ∫ ∇Nᵀ k ∇N dV (full rule), the capacity ∫ ρc N Nᵀ dV and a source ∫ N Q dV (mass rule), and the heat flux −k∇T at the integration points; heat::face_weights and face_film integrate surface flux and convection on tri3/tri6/quad4/quad8 faces. The solver (fieldcraft_solver::heat) numbers the temperatures with DofLayout::uniform(…, DofSet::TEMPERATURE) and assembles with the same sparse::add_block.

Thermal strain ε_th = α (T − T_ref) along the material axes enters a solid as the load Solid::thermal_load(d, &ThermalStrain { alpha, delta }) = ∫ Bᵀ D ε_th dV and its stress as Solid::stresses_with(d, u, Some(&thermal)) = D (ε − ε_th); delta is the nodal T − T_ref. 2D elements have the same pair (Plane::thermal_load, Plane::stresses_with): plane stress leaves εzz free (σzz = 0), plane strain holds it (σzz = −D₃·ε_th), axisymmetric elements strain the hoop direction too. Trusses (rod2, rod3) take the axial part Σ αₖ tₖ² ΔT along their tangent (Truss::thermal_load, Truss::stresses_with). Shells and beams will provide the same two functions for their strains; springs, dashpots and point masses have no thermal strain.

Explicit dynamics takes a step's temperature loads as a thermal strain rate per element, α ΔT / end time (the change reached linearly over the step, like a prescribed displacement, so there is no thermal shock), subtracted from each increment's strain increment before the constitutive update, so it feeds J2 and Johnson–Cook plasticity alike; isotropic expansion only (the rate is the same in any frame). Explicit trusses take no thermal strain yet (a warning says so).

Model side: Material::thermal (conductivity, specific_heat, expansion, each isotropic or per axis, and reference_temperature), BoundaryCondition::Temperature, loads SurfaceFlux, Convection, HeatSource (heat transfer) and Temperature (a temperature field in a static step), procedures HeatTransfer { transient } and CoupledTemperatureDisplacement. Commands: solve.heat, and solve.static for coupled steps. The deck writers refuse heat transfer by name for now.

Hybrid (mixed u–p) solids

A section with "formulation": {"hybrid": true} gives Formulation::Hybrid { pressures } (crate::hybrid): the isotropic elasticity splits as D = D_dev + K m mᵀ, the mean stress is an independent field discontinuous between elements (hex20: 1, x, y, z; hex8 and tet10: constant), and it is condensed out element by element, p̂ = K M⁻¹ Gᵀ u, leaving the stiffness K_dev + K G M⁻¹ Gᵀ. ElementDofs.internal.pressure counts the condensed unknowns; nothing global changes, and Cholesky still applies. Stresses are D_dev ε + p m with the element's pressure field. The choice of pressure space follows the inf-sup analysis (docs in crates/fem/src/hybrid.rs): Q2–P1disc and P2–P0 are stable; Q1–P0 does not lock but may show checkerboard pressures; linear tetrahedra and wedges lock even so and are refused. A thermal strain splits the same way: ∫ Bᵀ D_dev ε_th on the nodes and, through the pressure, K G M⁻¹ ∫ N_p tr ε_th, so a heated nearly incompressible solid takes its bulk stress −3KαΔT when confined. A continuous nodal pressure (DofKind::P, Taylor–Hood) needs a symmetric indefinite factorization and is planned.

Nonlinear static analysis (finite strain)

A step {"type": "nonlinear_static", …} (Procedure::NonlinearStatic, solved by solve.nonlinear, fieldcraft_solver::nonlinear) follows equilibrium in load increments with Newton–Raphson iterations. The step's load factor λ runs from 0 to 1: dead loads are λ times their values, a pressure follows the deformed faces with λ times its magnitude, and prescribed displacements are λ times theirs (a linear ramp).

Material point (fieldcraft_fem::constitutive). Given the deformation gradient F and the history committed at the end of the last converged increment, a material returns the Kirchhoff stress τ = J σ and the spatial tangent c of its Lie derivative (L_v τ = c : d), in Voigt order with engineering-shear columns, so c = D at F = I. This is the push-forward of the second Piola–Kirchhoff stress and its tangent: τ = F S Fᵀ, c_ijkl = F_iI F_jJ F_kK F_lL C_IJKL. History variables (MaterialModel::history_len per integration point) are read from the committed slice and written to a separate one, which the solver commits when the increment converges, so a cut-back increment leaves no trace. The materials:

Material in a nonlinear step	Model
Elastic (isotropic or orthotropic), no option	St Venant–Kirchhoff: S = D E, E = ½(FᵀF − I). Right for large rotations with small strains; it softens without bound in large compression.
Isotropic elastic with "hyperelastic": {"type": "neo_hookean"}	Compressible neo-Hookean with the isochoric–volumetric split, W = μ/2 (Ī₁ − 3) + κ/2 (J − 1)², μ and κ from E and ν (C₁₀ = μ/2, D₁ = 2/κ). Its small-strain limit is the linear material, which linear steps use.

"Isotropic elastic" keeps meaning the same thing everywhere: the neo-Hookean model is chosen only by the explicit option, and explicit dynamics refuses a hyperelastic material by name.

Elements (fieldcraft_fem::finite), total Lagrangian over the reference element with its own quadrature: f_a = ∫ Bₓ,aᵀ τ dV₀ and K_ab = ∫ (Bₓ,aᵀ c Bₓ,b + ∇ₓN_a·τ∇ₓN_b I) dV₀ with the current gradients ∇ₓN = F⁻ᵀ∇₀N (material plus geometric stiffness), identical term by term to ∫ B₀ᵀ S dV₀ and its linearisation. Every element's tangent is checked against central differences of its forces.

Element	Nonlinear steps
tet4, tet10, wedge6, hex8, hex20 (full), hex20 (reduced)	total Lagrangian
hybrid hex8, hex20, tet10	Simo–Taylor–Pister three-field form on the linear hybrid elements' pressure fields (θ the projection of J, p the projection of U′(θ)); needs the neo-Hookean option (St Venant–Kirchhoff has no volumetric split)
hex8 and quad4 incompatible modes	refused: their finite-strain form (enhanced assumed strain, Simo–Armero 1992) needs the modes as unknowns carried between iterations and hourglasses in large compression; use full integration, hybrid or quadratic elements
plane strain, axisymmetric (tri3, tri6, quad4, quad8)	total Lagrangian with F₃₃ = 1 or r/R (hoop column N/r, geometric term τ_θθ N_a N_b/r²)
plane stress	refused (σ_zz = 0 needs a local iteration for the thickness)
trusses, shells, beams	refused (large rotations of structural elements come later in M15)

Without nonlinear geometry ("nonlinear_geometry": false) the elements are the linear ones: every material is elastic today, so such a step is the linear solution reached by Newton (the path the next materials, plasticity first, take with small strains).

Follower pressure (fieldcraft_fem::face::follower_pressure, follower_edge_pressure): fᵢ = −p ∫ Nᵢ (x,ξ × x,η) dξ dη on the current face, with its load stiffness K = p ∫ N_a (N_b,η [x,ξ]× − N_b,ξ [x,η]×). The load stiffness is not symmetric. The solver keeps its symmetric part and leaves the skew part out, so the tangent stays symmetric for the supernodal Cholesky factorization: the residual is exact (the converged answer does not change), and over a closed surface or one whose boundary is held the skew parts of neighbouring faces cancel (Schweizerhof & Ramm 1984), so nothing is lost there; elsewhere Newton may converge linearly instead of quadratically, and the result says how large the left-out part is. An LU factorization would double the factor's memory and lose the symmetric solver for a gain only in the rate of convergence of open loaded surfaces; it is revisited with the symmetric indefinite (LDLᵀ) factorization that arc-length and Taylor–Hood need. Without nonlinear geometry a pressure is a dead load.

Newton–Raphson with the consistent tangent refactored every iteration, an optional line search (off by default: on large elastic rotations the full step converges in fewer iterations), and convergence when the force residual (∞-norm over the largest force seen), the displacement correction (over the increment's largest displacement) and the energy (over the first iteration's) are all below their tolerances, or the residual is round-off. Each increment starts from the linear extrapolation of the last one. Automatic incrementation grows an increment by half after two fast ones (≤ 4 iterations), shrinks it after a slow one (> 8), and cuts it to a quarter when an element inverts, the residual is not finite or grows twice, the tangent is not positive definite, or the iterations run out; the step stops with the load factor and the reason at the minimum increment, the increment cap, or Limits::max_newton_iterations.

Results: the final displacements and reactions, the Cauchy stress and the logarithmic strain (ln V, engineering shears) at the integration points and averaged at the nodes, every increment's load factor, iterations and residuals, and a displacement frame per converged increment (results.frame, the slider in the results panel).

Couplings stay the linear constraints of the linear solver (small rotations of the reference points), with a warning; springs act along their fixed directions. Temperature and centrifugal loads are refused by name in nonlinear steps.

Decks: Abaqus *STEP, NLGEOM=YES|NO, INC= with *STATIC (initial, period 1, minimum, maximum) and *CONTROLS, PARAMETERS=LINE SEARCH; the neo-Hookean option as *HYPERELASTIC, NEO HOOKE. LS-DYNA: the ramped implicit run with full Newton (ilimit 1, maxref) and *CONTROL_IMPLICIT_AUTO; hyperelastic materials in nonlinear steps and steps without nonlinear geometry are refused. Nastran refuses nonlinear steps (SOL 400) by name. Iteration caps and tolerances are left to each solver's own controls, whose norms differ from ours.

Verification (crates/solver/tests/nonlinear.rs, crates/fem/src/{constitutive,finite,face}.rs, and tests/oracle/nl-*): consistent tangents of both materials and of every element (central differences); uniaxial stretch to λ = 2 of a neo-Hookean block in every element kind, plane strain and axisymmetric (closed form); St Venant–Kirchhoff uniaxial stretch; the nonlinear patch test on a distorted mesh (all solids, hybrid ones included); a 90° rigid rotation (zero stress); the elastica (Bisshopp–Drucker) at P L²/(E I) = 1, 2, 5, 10 with hex20 (0.2–0.5 %); the inflation of a thick neo-Hookean tube by a follower pressure (Rivlin's closed form); Newton's quadratic convergence; the linear limit; a finite-strain Cook's membrane is pending its transcribed reference; and the failures (an element inverted by crushing stops at the minimum increment, naming it; the caps; what is refused).