FEM foundation B: heat transfer and thermal stress (solids, 2D, explicit), coupled temperature–displacement, solve.heat (#45)

* Heat transfer and thermal stress: thermal materials, temperature BCs, heat loads, heat transfer and coupled temperature–displacement steps, solve.heat

Model: Material.thermal (conductivity, specific heat, expansion, isotropic or
per axis; reference temperature), BoundaryCondition::Temperature (dof 11),
loads SurfaceFlux, Convection, HeatSource and Temperature (a field in static
steps), procedures HeatTransfer { transient (θ-method) } and
CoupledTemperatureDisplacement; checks for missing conductivity/specific
heat/expansion. fem: heat conduction kernels (conductivity, capacity,
source, face flux/film, heat flux) on the DofLayout (one T per node), thermal
strain load and stress D(ε − ε_th) for every solid formulation incl.
incompatible modes. Solver: steady and transient heat transfer, thermal loads
in static steps, coupled steps (exact one-way coupling). Engine: solve.heat,
temperature/heat-flux result fields, probes. Decks refuse heat transfer by
name. Analytic verification: slab, source+convection/flux, thick cylinder,
erfc semi-infinite solid, free/constrained expansion, bimetal strip, coupled
step.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* thermal: strains in 2D elements and explicit dynamics

Plane::thermal_load and stresses_with for plane stress (εzz free), plane
strain (σzz = −D₃·ε_th) and axisymmetric elements, incompatible modes
included; the solver's 2D kernels take them. Explicit dynamics turns a
step's temperature loads into a per-element thermal strain rate (reached
linearly over the step) subtracted from each strain increment, so it feeds
J2 and Johnson–Cook updates. Verified: free expansion in every 2D theory
and element, a biaxially held plate, an explicit free bar (ΔL = αΔT L,
stress-free) and a clamped bar at −EαΔT, capped at the yield stress with
perfect plasticity.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* clippy

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* thermal: trusses take the axial thermal strain; keep the roadmap's M7 line

Truss::thermal_load and stresses_with project the thermal strain on the
member's tangent (Σ αₖ tₖ² ΔT); the static solver applies them to rod2/rod3
(#38). Explicit trusses take none yet and a warning says so.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
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 <noreply@anthropic.com>

crates/fem/src/element.rs

@@ -195,6 +195,26 @@ struct Gradients {
     det: f64,
 }
 
+/// A thermal strain field over an element: the expansion coefficients along X, Y, Z and the
+/// nodal temperature changes T − T_ref.
+#[derive(Debug, Clone, PartialEq)]
+pub struct ThermalStrain {
+    pub alpha: [f64; 3],
+    pub delta: Vec<f64>,
+}
+
+impl ThermalStrain {
+    pub(crate) fn check(&self, nodes: usize) -> Result<()> {
+        if self.delta.len() != nodes {
+            return Err(Error::Invalid(format!("{} temperatures for a {nodes}-node element", self.delta.len())));
+        }
+        if self.alpha.iter().chain(&self.delta).any(|v| !v.is_finite()) {
+            return Err(Error::Invalid("temperatures and expansion coefficients must be finite".into()));
+        }
+        Ok(())
+    }
+}
+
 /// A solid element: its type, its node coordinates and its formulation.
 #[derive(Debug, Clone, Copy)]
 pub struct Solid<'a> {
@@ -423,16 +443,83 @@ impl<'a> Solid<'a> {
         Ok(f)
     }
 
+    /// The thermal strain α ΔT at a point with shape values `n` (xx, yy, zz; no shear).
+    fn thermal_strain(thermal: &ThermalStrain, n: &[f64]) -> [f64; 6] {
+        let dt: f64 = n.iter().zip(&thermal.delta).map(|(a, b)| a * b).sum();
+        [thermal.alpha[0] * dt, thermal.alpha[1] * dt, thermal.alpha[2] * dt, 0.0, 0.0, 0.0]
+    }
+
+    /// The equivalent nodal forces of a thermal strain, ∫ Bᵀ D ε_th dV, with the stiffness rule
+    /// (the initial-strain load of Zienkiewicz, Taylor & Zhu 2013, §6.2.3); incompatible modes
+    /// are condensed out as in the stiffness, f = f_u − K_ua K_aa⁻¹ f_a.
+    pub fn thermal_load(&self, d: &Elasticity, thermal: &ThermalStrain) -> Result<Vec<[f64; 3]>> {
+        thermal.check(self.ty.nodes)?;
+        if self.formulation == Formulation::IncompatibleModes {
+            let IncompatibleParts { kua, kaa, j0, .. } = self.incompatible_parts(d, false)?;
+            let (mut f, mut fa) = self.thermal_forces(d, thermal, Some(&j0))?;
+            Cholesky::new(&kaa)?.solve(&mut fa);
+            for (i, fi) in f.iter_mut().flatten().enumerate() {
+                *fi -= (0..INCOMPATIBLE_DOFS).map(|a| kua.get(i * INCOMPATIBLE_DOFS + a).copied().unwrap_or(0.0) * fa[a]).sum::<f64>();
+            }
+            return Ok(f);
+        }
+        Ok(self.thermal_forces(d, thermal, None)?.0)
+    }
+
+    /// ∫ Bᵀ D ε_th dV on the nodes and, with the centre Jacobian `j0`, on the incompatible modes.
+    fn thermal_forces(
+        &self,
+        d: &Elasticity,
+        thermal: &ThermalStrain,
+        j0: Option<&([[f64; 3]; 3], f64)>,
+    ) -> Result<(Vec<[f64; 3]>, [f64; INCOMPATIBLE_DOFS])> {
+        let mut f = vec![[0.0; 3]; self.ty.nodes];
+        let mut fa = [0.0; INCOMPATIBLE_DOFS];
+        let add = |out: &mut [f64], g: &[f64; 3], sig: &[f64; 6], w: f64| {
+            let b = b_block(g);
+            for (c, o) in out.iter_mut().enumerate() {
+                *o += w * (0..6).map(|s| b[s][c] * sig[s]).sum::<f64>();
+            }
+        };
+        for p in self.rule().points() {
+            let g = self.gradients(p.xi)?;
+            let w = g.det * p.weight;
+            let sig = stress(d, &Self::thermal_strain(thermal, g.shape.values()));
+            for (fi, gi) in f.iter_mut().zip(g.dx.get(..self.ty.nodes).unwrap_or_default()) {
+                add(fi, gi, &sig, w);
+            }
+            if let Some(j0) = j0 {
+                for (k, gk) in self.incompatible_gradients(p.xi, g.det, j0).iter().enumerate() {
+                    add(fa.get_mut(3 * k..3 * k + 3).unwrap_or_default(), gk, &sig, w);
+                }
+            }
+        }
+        Ok((f, fa))
+    }
+
     /// Strain and stress at every integration point for the nodal displacements `u`.
     pub fn stresses(&self, d: &Elasticity, u: &[[f64; 3]]) -> Result<Vec<PointStress>> {
+        self.stresses_with(d, u, None)
+    }
+
+    /// Strain and stress at every integration point for the nodal displacements `u` and, if
+    /// given, a thermal strain: the stress is D (ε − ε_th), the strain reported the total ε.
+    pub fn stresses_with(&self, d: &Elasticity, u: &[[f64; 3]], thermal: Option<&ThermalStrain>) -> Result<Vec<PointStress>> {
         if u.len() != self.ty.nodes {
             return Err(Error::Invalid(format!("{} displacements for a {}-node element", u.len(), self.ty.nodes)));
         }
-        // Incompatible modes: recover the internal amplitudes α = −K_aa⁻¹ K_au u.
+        if let Some(t) = thermal {
+            t.check(self.ty.nodes)?;
+        }
+        // Incompatible modes: recover the internal amplitudes from K_aa α + K_au u = f_a (the
+        // thermal load on the modes; zero without one).
         let alpha = if self.formulation == Formulation::IncompatibleModes {
             let IncompatibleParts { kua, kaa, j0, .. } = self.incompatible_parts(d, false)?;
             let flat: Vec<f64> = u.iter().flatten().copied().collect();
-            let mut rhs = [0.0; INCOMPATIBLE_DOFS];
+            let mut rhs = match thermal {
+                Some(t) => self.thermal_forces(d, t, Some(&j0))?.1,
+                None => [0.0; INCOMPATIBLE_DOFS],
+            };
             for (i, ui) in flat.iter().enumerate() {
                 for (a, r) in rhs.iter_mut().enumerate() {
                     *r -= kua.get(i * INCOMPATIBLE_DOFS + a).copied().unwrap_or(0.0) * ui;
@@ -471,7 +558,14 @@ impl<'a> Solid<'a> {
                     xyz[k] += n * x[k];
                 }
             }
-            out.push(PointStress { xyz, strain, stress: stress(d, &strain) });
+            let sigma = match thermal {
+                Some(t) => {
+                    let eth = Self::thermal_strain(t, g.shape.values());
+                    stress(d, &[0, 1, 2, 3, 4, 5].map(|k| strain[k] - eth[k]))
+                }
+                None => stress(d, &strain),
+            };
+            out.push(PointStress { xyz, strain, stress: sigma });
         }
         Ok(out)
     }

crates/fem/src/heat.rs

@@ -0,0 +1,233 @@
+//! Heat conduction kernels: conductivity, heat capacity, internal heat sources, surface fluxes and
+//! convection, and the heat flux at the integration points.
+//!
+//! The Galerkin form of transient conduction ρc Ṫ = ∇·(k∇T) + Q with flux and convection on the
+//! boundary (Zienkiewicz, Taylor & Zhu 2013, §3.2 and ch. 7 of vol. 1; Hughes 2000, ch. 1, §8.1)
+//! gives, element by element, with one temperature per node ([`ElementDofs`] `TEMPERATURE`):
+//! - conductivity K = ∫ ∇Nᵀ k ∇N dV, k the diagonal conductivity in the global axes (the
+//!   material axes until orientations land), integrated with the element's full rule;
+//! - capacity C = ∫ ρc N Nᵀ dV with the mass rule (consistent);
+//! - source f = ∫ N Q dV;
+//! - on faces: a flux q into the body gives f = ∫ N q dA; convection with film coefficient h to a
+//!   sink at T∞ adds h ∫ N Nᵀ dA to K and h T∞ ∫ N dA to f.
+//!
+//! Faces are integrated on their own (possibly curved) geometry with the model's face rule
+//! (the one pressures use), exact for N Nᵀ on flat faces.
+
+use modelcraft_model::ElementKind;
+
+use crate::basis::Basis;
+use crate::dense::Matrix;
+use crate::dof::{DofSet, ElementDofs};
+use crate::element::{ElementType, element_type};
+use crate::quadrature::Rule;
+use crate::{Error, Result};
+
+/// What a heat conduction element carries: one temperature per node.
+pub const HEAT_DOFS: ElementDofs = ElementDofs::new(DofSet::TEMPERATURE);
+
+/// The heat flux −k∇T at one integration point.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct PointFlux {
+    pub xyz: [f64; 3],
+    pub flux: [f64; 3],
+}
+
+/// A solid element conducting heat.
+#[derive(Debug, Clone, Copy)]
+pub struct Conduction<'a> {
+    pub ty: &'static ElementType,
+    basis: Basis,
+    xyz: &'a [[f64; 3]],
+}
+
+impl<'a> Conduction<'a> {
+    pub fn new(kind: ElementKind, xyz: &'a [[f64; 3]]) -> Result<Conduction<'a>> {
+        let ty = element_type(kind)?;
+        if xyz.len() != ty.nodes {
+            return Err(Error::Invalid(format!("a {} element has {} nodes, got {}", kind.name(), ty.nodes, xyz.len())));
+        }
+        modelcraft_model::check_coordinates(xyz.as_flattened()).map_err(|e| Error::Invalid(e.to_string()))?;
+        Ok(Conduction { ty, basis: Basis::of(kind), xyz })
+    }
+
+    /// Visit the points of `rule`: shape values, physical gradients and the weight w det J.
+    fn each(&self, rule: Rule, mut visit: impl FnMut(&[f64], &[[f64; 3]], f64)) -> Result<()> {
+        for p in rule.points() {
+            let s = self.basis.eval(p.xi);
+            let g = self.basis.gradients(&s, self.xyz)?;
+            visit(s.values(), g.gradients(), g.det * p.weight);
+        }
+        Ok(())
+    }
+
+    /// K = ∫ ∇Nᵀ k ∇N dV, n × n.
+    pub fn conductivity(&self, k: [f64; 3]) -> Result<Matrix> {
+        let mut m = Matrix::zeros(self.ty.nodes);
+        self.each(self.ty.full, |_, g, w| {
+            for (i, gi) in g.iter().enumerate() {
+                for (j, gj) in g.iter().enumerate() {
+                    m.add(i, j, w * (k[0] * gi[0] * gj[0] + k[1] * gi[1] * gj[1] + k[2] * gi[2] * gj[2]));
+                }
+            }
+        })?;
+        Ok(m)
+    }
+
+    /// C = ∫ ρc N Nᵀ dV, n × n; `capacity` is ρc (energy per volume per temperature).
+    pub fn capacity(&self, capacity: f64) -> Result<Matrix> {
+        let mut m = Matrix::zeros(self.ty.nodes);
+        self.each(self.ty.mass, |n, _, w| {
+            for (i, a) in n.iter().enumerate() {
+                for (j, b) in n.iter().enumerate() {
+                    m.add(i, j, w * capacity * a * b);
+                }
+            }
+        })?;
+        Ok(m)
+    }
+
+    /// f = ∫ N Q dV for a uniform heat source Q per volume.
+    pub fn source(&self, q: f64) -> Result<Vec<f64>> {
+        let mut f = vec![0.0; self.ty.nodes];
+        self.each(self.ty.mass, |n, _, w| {
+            for (fi, ni) in f.iter_mut().zip(n) {
+                *fi += w * q * ni;
+            }
+        })?;
+        Ok(f)
+    }
+
+    /// The heat flux −k∇T at the full rule's points for nodal temperatures `t`.
+    pub fn flux(&self, k: [f64; 3], t: &[f64]) -> Result<Vec<PointFlux>> {
+        if t.len() != self.ty.nodes {
+            return Err(Error::Invalid(format!("{} temperatures for a {}-node element", t.len(), self.ty.nodes)));
+        }
+        let mut out = Vec::new();
+        self.each(self.ty.full, |n, g, _| {
+            let (mut grad, mut xyz) = ([0.0; 3], [0.0; 3]);
+            for ((gi, ti), (ni, xi)) in g.iter().zip(t).zip(n.iter().zip(self.xyz)) {
+                for c in 0..3 {
+                    grad[c] += gi[c] * ti;
+                    xyz[c] += ni * xi[c];
+                }
+            }
+            out.push(PointFlux { xyz, flux: [0, 1, 2].map(|c| -k[c] * grad[c]) });
+        })?;
+        Ok(out)
+    }
+}
+
+/// Visit the points of a face (node coordinates in the face's node order): shape values and
+/// the area weight, with the model's face rule ([`modelcraft_model::nodal::integrate_face`], the
+/// one pressures and tractions use).
+fn each_face_point(xyz: &[[f64; 3]], mut visit: impl FnMut(&[f64], f64)) -> Result<()> {
+    let mut collapsed = false;
+    modelcraft_model::nodal::integrate_face(xyz, |n, nda| {
+        let area = (nda[0] * nda[0] + nda[1] * nda[1] + nda[2] * nda[2]).sqrt();
+        collapsed |= !(area > 0.0 && area.is_finite());
+        visit(n, area);
+    })
+    .map_err(|e| Error::Invalid(e.to_string()))?;
+    if collapsed {
+        return Err(Error::Degenerate("a face is collapsed (zero area)".into()));
+    }
+    Ok(())
+}
+
+/// ∫ N dA for every node of a face: the nodal shares of a uniform flux.
+pub fn face_weights(xyz: &[[f64; 3]]) -> Result<Vec<f64>> {
+    let mut f = vec![0.0; xyz.len()];
+    each_face_point(xyz, |n, w| {
+        for (fi, ni) in f.iter_mut().zip(n) {
+            *fi += w * ni;
+        }
+    })?;
+    Ok(f)
+}
+
+/// h ∫ N Nᵀ dA for a face: the conductivity convection adds.
+pub fn face_film(xyz: &[[f64; 3]], film: f64) -> Result<Matrix> {
+    let mut m = Matrix::zeros(xyz.len());
+    each_face_point(xyz, |n, w| {
+        for (i, a) in n.iter().enumerate() {
+            for (j, b) in n.iter().enumerate() {
+                m.add(i, j, w * film * a * b);
+            }
+        }
+    })?;
+    Ok(m)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::shape::{HEX_CORNERS, TET_CORNERS};
+
+    fn box_hex(a: [f64; 3]) -> Vec<[f64; 3]> {
+        HEX_CORNERS.iter().map(|c| [0, 1, 2].map(|k| (c[k] + 1.0) / 2.0 * a[k])).collect()
+    }
+
+    #[test]
+    fn conductivity_rows_sum_to_zero_and_reproduce_a_linear_field() {
+        // A uniform gradient G in a 2 × 1 × 3 box: K T = boundary flux, zero at interior
+        // balance; here check K·1 = 0 and Tᵀ K T = k |G|² V.
+        let xyz = box_hex([2.0, 1.0, 3.0]);
+        for kind in [ElementKind::Hex8] {
+            let e = Conduction::new(kind, &xyz).unwrap();
+            let k = [2.0, 3.0, 5.0];
+            let m = e.conductivity(k).unwrap();
+            let ones = vec![1.0; 8];
+            assert!(m.mul_vec(&ones).iter().all(|v| v.abs() < 1e-12));
+            let g = [0.3, -0.7, 0.2];
+            let t: Vec<f64> = xyz.iter().map(|x| g[0] * x[0] + g[1] * x[1] + g[2] * x[2]).collect();
+            let kt = m.mul_vec(&t);
+            let energy: f64 = kt.iter().zip(&t).map(|(a, b)| a * b).sum();
+            let exact = (k[0] * g[0] * g[0] + k[1] * g[1] * g[1] + k[2] * g[2] * g[2]) * 6.0;
+            assert!((energy - exact).abs() < 1e-12 * exact, "{energy} vs {exact}");
+            let fl = e.flux(k, &t).unwrap();
+            for p in fl {
+                for c in 0..3 {
+                    assert!((p.flux[c] + k[c] * g[c]).abs() < 1e-13);
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn capacity_and_source_integrate_the_volume() {
+        let tet: Vec<[f64; 3]> = TET_CORNERS.iter().map(|c| [2.0 * c[0], 3.0 * c[1], c[2]]).collect();
+        let e = Conduction::new(ElementKind::Tet4, &tet).unwrap();
+        let v = 1.0;
+        let c = e.capacity(4.0).unwrap();
+        assert!((c.data.iter().sum::<f64>() - 4.0 * v).abs() < 1e-13);
+        assert!((c.get(0, 0) - 4.0 * v / 10.0).abs() < 1e-13 && (c.get(0, 1) - 4.0 * v / 20.0).abs() < 1e-13);
+        let f = e.source(3.0).unwrap();
+        assert!(f.iter().all(|x| (x - 0.75).abs() < 1e-13), "{f:?}");
+        assert!(Conduction::new(ElementKind::Tet4, &tet[..3]).is_err());
+        assert!(e.flux([1.0; 3], &[0.0; 3]).is_err());
+    }
+
+    #[test]
+    fn face_weights_and_film() {
+        let quad = [[0., 0., 0.], [2., 0., 0.], [2., 3., 0.], [0., 3., 0.]];
+        let w = face_weights(&quad).unwrap();
+        assert!(w.iter().all(|x| (x - 1.5).abs() < 1e-14));
+        let m = face_film(&quad, 2.0).unwrap();
+        // h A / 9 on the diagonal of a bilinear rectangle, sum h A.
+        assert!((m.data.iter().sum::<f64>() - 12.0).abs() < 1e-12);
+        assert!((m.get(0, 0) - 2.0 * 6.0 / 9.0).abs() < 1e-13);
+        let mut t6 = vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]];
+        for (a, b) in [(0, 1), (1, 2), (2, 0)] {
+            t6.push([0, 1, 2].map(|k| (t6[a][k] + t6[b][k]) / 2.0));
+        }
+        let w = face_weights(&t6).unwrap();
+        for (i, x) in w.iter().enumerate() {
+            let want = if i < 3 { 0.0 } else { 0.5 / 3.0 };
+            assert!((x - want).abs() < 1e-14, "{i}: {x}");
+        }
+        assert!(face_weights(&quad[..2]).is_err());
+        assert!(face_weights(&[[0.0; 3]; 4]).is_err());
+        assert!(face_film(&[[f64::NAN, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], 1.0).is_err());
+    }
+}

crates/fem/src/lib.rs

@@ -22,6 +22,8 @@
 //! - [`beam`]: two-node Timoshenko beams with profile section properties;
 //! - [`frame`]: local axes of shells and beams and the rotation of their matrices;
 //! - [`face`]: consistent nodal forces of pressure on element faces;
+//! - [`heat`]: heat conduction (conductivity, capacity, sources, surface flux and convection,
+//!   heat flux); thermal strain loads and stresses are in [`element`];
 //! - [`uniform`]: uniform-strain (one-point) hexahedron and tetrahedron kinematics, hourglass
 //!   shape vectors and element frequency bounds for explicit dynamics;
 //! - [`rotation`]: incrementally objective stress updates (Hughes–Winget Jaumann, Flanagan–Taylor
@@ -47,6 +49,7 @@ pub mod dof;
 pub mod element;
 pub mod face;
 pub mod frame;
+pub mod heat;
 pub mod jet;
 pub mod material;
 pub mod plane;
@@ -60,7 +63,7 @@ pub mod uniform;
 
 pub use basis::Basis;
 pub use dof::{DofKind, DofLayout, DofSet, ElementDofs};
-pub use element::{ELEMENT_TYPES, ElementType, Family, Formulation, PointStress, Solid, element_type};
+pub use element::{ELEMENT_TYPES, ElementType, Family, Formulation, PointStress, Solid, ThermalStrain, element_type};
 pub use material::{Elasticity, elasticity, principal, von_mises};
 
 #[derive(Debug, Clone, PartialEq, thiserror::Error)]

crates/fem/src/plane.rs

@@ -44,7 +44,7 @@ use modelcraft_model::{ElementKind, Extent, SectionKind};
 use crate::basis::Basis;
 use crate::dense::{Cholesky, Matrix};
 use crate::dof::{DofSet, ElementDofs, Internal};
-use crate::element::{ElementType, Family, Formulation, PointStress, extrapolation, lump, symmetrize};
+use crate::element::{ElementType, Family, Formulation, PointStress, ThermalStrain, extrapolation, lump, symmetrize};
 use crate::material::Elasticity;
 use crate::quadrature::{Point, Rule};
 use crate::shape::{MAX_NODES, Shape};
@@ -471,20 +471,98 @@ impl<'a> Plane<'a> {
         Ok(f)
     }
 
+    /// The thermal strain α ΔT at a point with shape values `n` (xx, yy, zz, xy; zz is the hoop
+    /// direction of an axisymmetric element).
+    fn thermal_strain(thermal: &ThermalStrain, n: &[f64]) -> [f64; 4] {
+        let dt: f64 = n.iter().zip(&thermal.delta).map(|(a, b)| a * b).sum();
+        [thermal.alpha[0] * dt, thermal.alpha[1] * dt, thermal.alpha[2] * dt, 0.0]
+    }
+
+    /// ∫ Bᵀ D ε_th dV on the nodes' (u, v) and on the incompatible modes.
+    fn thermal_forces(
+        d2: &Elasticity2,
+        thermal: &ThermalStrain,
+        points: &[(Point, At)],
+        modes: Option<&Vec<[Column; 2]>>,
+        nodes: usize,
+    ) -> (Vec<f64>, Vec<f64>) {
+        let mut fu = vec![0.0; 2 * nodes];
+        let mut fa = vec![0.0; INCOMPATIBLE_DOFS];
+        for (i, (_, a)) in points.iter().enumerate() {
+            let eps = Self::thermal_strain(thermal, a.shape.values());
+            let s: [f64; 4] = std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * eps[c]).sum());
+            // Bᵀ σ for one strain column g: u takes εxx, εzz and γxy; v takes εyy and γxy.
+            let add = |out: &mut [f64], k: usize, g: &Column| {
+                if let Some(x) = out.get_mut(2 * k) {
+                    *x += a.dv * (g[0] * s[0] + g[2] * s[2] + g[1] * s[3]);
+                }
+                if let Some(x) = out.get_mut(2 * k + 1) {
+                    *x += a.dv * (g[1] * s[1] + g[0] * s[3]);
+                }
+            };
+            for (k, g) in a.g.iter().take(nodes).enumerate() {
+                add(&mut fu, k, g);
+            }
+            if let Some(cols) = modes.and_then(|m| m.get(i)) {
+                for (k, g) in cols.iter().enumerate() {
+                    add(&mut fa, k, g);
+                }
+            }
+        }
+        (fu, fa)
+    }
+
+    /// The equivalent nodal forces (u, v) of a thermal strain, ∫ Bᵀ D ε_th dV with the
+    /// theory's elasticity (plane stress leaves εzz free, plane strain holds it, axisymmetric
+    /// strains the hoop direction too); incompatible modes are condensed as in the stiffness.
+    pub fn thermal_load(&self, d: &Elasticity, thermal: &ThermalStrain) -> Result<Vec<[f64; 2]>> {
+        thermal.check(self.ty.nodes)?;
+        let d2 = self.theory.elasticity(d)?;
+        let n = 2 * self.ty.nodes;
+        let mut fu = if self.formulation == Formulation::IncompatibleModes {
+            let (full, points, modes) = self.full_stiffness(&d2)?;
+            let (mut fu, mut fa) = Self::thermal_forces(&d2, thermal, &points, modes.as_ref(), self.ty.nodes);
+            Cholesky::new(&sub_block(&full, n, n, INCOMPATIBLE_DOFS))?.solve(&mut fa);
+            for (i, x) in fu.iter_mut().enumerate() {
+                *x -= (0..INCOMPATIBLE_DOFS).map(|a| full.get(i, n + a) * fa.get(a).copied().unwrap_or(0.0)).sum::<f64>();
+            }
+            fu
+        } else {
+            Self::thermal_forces(&d2, thermal, &self.points()?, None, self.ty.nodes).0
+        };
+        fu.resize(n, 0.0);
+        Ok(fu.as_chunks::<2>().0.to_vec())
+    }
+
     /// Strain and stress at every integration point for the nodal displacements `u` (u, v).
     pub fn stresses(&self, d: &Elasticity, u: &[[f64; 2]]) -> Result<Vec<PointStress>> {
+        self.stresses_with(d, u, None)
+    }
+
+    /// Strain and stress at every integration point for the nodal displacements `u` (u, v) and,
+    /// if given, a thermal strain: the stress is D (ε − ε_th), the strain reported the total.
+    pub fn stresses_with(&self, d: &Elasticity, u: &[[f64; 2]], thermal: Option<&ThermalStrain>) -> Result<Vec<PointStress>> {
         if u.len() != self.ty.nodes {
             return Err(Error::Invalid(format!("{} displacements for a {}-node element", u.len(), self.ty.nodes)));
         }
+        if let Some(t) = thermal {
+            t.check(self.ty.nodes)?;
+        }
         let d2 = self.theory.elasticity(d)?;
         let (full, points, modes) =
             if self.formulation == Formulation::IncompatibleModes { self.full_stiffness(&d2)? } else { (Matrix::zeros(0), self.points()?, None) };
-        // Incompatible modes: α = −K_aa⁻¹ K_au u.
+        // Incompatible modes: K_aa α = f_a − K_au u (f_a the thermal load on the modes).
         let alpha = if modes.is_some() {
             let n = 2 * self.ty.nodes;
             let flat: Vec<f64> = u.iter().flatten().copied().collect();
             let mut rhs: Vec<f64> =
                 (0..INCOMPATIBLE_DOFS).map(|a| -flat.iter().enumerate().map(|(i, x)| full.get(n + a, i) * x).sum::<f64>()).collect();
+            if let Some(t) = thermal {
+                let (_, fa) = Self::thermal_forces(&d2, t, &points, modes.as_ref(), self.ty.nodes);
+                for (r, f) in rhs.iter_mut().zip(fa) {
+                    *r += f;
+                }
+            }
             Cholesky::new(&sub_block(&full, n, n, INCOMPATIBLE_DOFS))?.solve(&mut rhs);
             rhs
         } else {
@@ -507,11 +585,19 @@ impl<'a> Plane<'a> {
                     add(g, [alpha.get(2 * k).copied().unwrap_or(0.0), alpha.get(2 * k + 1).copied().unwrap_or(0.0)]);
                 }
             }
+            let eth = thermal.map(|t| Self::thermal_strain(t, a.shape.values()));
             if let Theory::PlaneStress { .. } = self.theory {
-                e[2] = -(d[2][0] * e[0] + d[2][1] * e[1] + d[2][3] * e[3]) / d[2][2];
+                e[2] = match eth {
+                    // σzz = 0 with the thermal strain: εzz = εth,zz − D₃·(ε − εth)/D₃₃ in plane.
+                    Some(t) => t[2] - (d[2][0] * (e[0] - t[0]) + d[2][1] * (e[1] - t[1]) + d[2][3] * e[3]) / d[2][2],
+                    None => -(d[2][0] * e[0] + d[2][1] * e[1] + d[2][3] * e[3]) / d[2][2],
+                };
             }
             // Plane stress: D̃ has no zz row or column, so the εzz just recovered adds nothing.
-            let s: [f64; 4] = std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * e[c]).sum());
+            let s: [f64; 4] = match eth {
+                Some(t) => std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * (e[c] - t[c])).sum()),
+                None => std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * e[c]).sum()),
+            };
             out.push(PointStress { xyz: a.x, strain: [e[0], e[1], e[2], e[3], 0.0, 0.0], stress: [s[0], s[1], s[2], s[3], 0.0, 0.0] });
         }
         Ok(out)

crates/fem/src/tests.rs

@@ -208,3 +208,33 @@ fn section_formulations() {
     assert!(Formulation::from_section(ElementKind::Tet4, &hybrid).is_err());
     assert!(Formulation::from_section(ElementKind::Quad4, &full).is_err());
 }
+
+#[test]
+fn uniform_thermal_strain_is_free_expansion_and_its_load_balances_it() {
+    // A uniform ΔT with orthotropic α: the free expansion u = diag(α ΔT) x is stress-free, and
+    // the thermal load equals K u for it (the element is in equilibrium with no other load).
+    let d = orthotropic([1000.0, 800.0, 600.0], [0.3, 0.25, 0.2], [400.0, 300.0, 350.0]).unwrap();
+    let alpha = [1.2e-5, 2.0e-5, 0.7e-5];
+    for (kind, f) in CASES {
+        let xyz = distorted(kind, true);
+        let solid = Solid::new(kind, &xyz, f).unwrap();
+        let th = crate::ThermalStrain { alpha, delta: vec![80.0; xyz.len()] };
+        let u: Vec<[f64; 3]> = xyz.iter().map(|x| [0, 1, 2].map(|k| alpha[k] * 80.0 * x[k])).collect();
+        for p in solid.stresses_with(&d, &u, Some(&th)).unwrap() {
+            assert!(p.stress.iter().all(|s| s.abs() < 1e-9), "{kind:?} {f:?}: {:?}", p.stress);
+            for k in 0..3 {
+                assert!((p.strain[k] - alpha[k] * 80.0).abs() < 1e-12);
+            }
+        }
+        let load = solid.thermal_load(&d, &th).unwrap();
+        let ku = solid.stiffness(&d).unwrap().mul_vec(&u.iter().flatten().copied().collect::<Vec<_>>());
+        let scale = ku.iter().fold(0.0f64, |m, v| m.max(v.abs()));
+        for (a, b) in load.iter().flatten().zip(&ku) {
+            assert!((a - b).abs() < 1e-9 * scale, "{kind:?} {f:?}: {a} vs {b}");
+        }
+        // Without a thermal strain the two stress paths agree bit for bit.
+        assert_eq!(solid.stresses_with(&d, &u, None).unwrap(), solid.stresses(&d, &u).unwrap());
+        assert!(solid.thermal_load(&d, &crate::ThermalStrain { alpha, delta: vec![1.0; 2] }).is_err());
+        assert!(solid.thermal_load(&d, &crate::ThermalStrain { alpha: [f64::NAN; 3], delta: vec![1.0; xyz.len()] }).is_err());
+    }
+}

crates/fem/src/truss.rs

@@ -29,7 +29,7 @@ use modelcraft_model::{Elastic, ElementKind};
 use crate::basis::Basis;
 use crate::dense::Matrix;
 use crate::dof::{DofSet, ElementDofs};
-use crate::element::PointStress;
+use crate::element::{PointStress, ThermalStrain};
 use crate::quadrature::gauss_legendre;
 use crate::{Error, Result};
 
@@ -196,12 +196,51 @@ impl<'a> Truss<'a> {
 
     /// Strain and stress at the integration points as global tensors (see the module docs).
     pub fn stresses(&self, young: f64, u: &[[f64; 3]]) -> Result<Vec<PointStress>> {
-        let strains = self.axial_strains(u)?;
+        self.stresses_with(young, u, None)
+    }
+
+    /// The axial thermal strain at each integration point: the thermal strain tensor
+    /// diag(α) ΔT along the tangent, Σ αₖ tₖ² ΔT, with ΔT interpolated from the nodes.
+    fn thermal_strains(&self, thermal: &ThermalStrain) -> Result<Vec<f64>> {
+        thermal.check(self.nodes())?;
         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) })
+            .map(|p| {
+                let dt: f64 = p.n.iter().zip(&thermal.delta).map(|(a, b)| a * b).sum();
+                (0..3).map(|k| thermal.alpha[k] * p.tangent[k] * p.tangent[k]).sum::<f64>() * dt
+            })
+            .collect())
+    }
+
+    /// The equivalent nodal forces of a thermal strain, ∫ Bᵀ E A ε_th ds.
+    pub fn thermal_load(&self, young: f64, area: f64, thermal: &ThermalStrain) -> Result<Vec<[f64; 3]>> {
+        let strains = self.thermal_strains(thermal)?;
+        let mut f = vec![[0.0; 3]; self.nodes()];
+        for (p, e) in self.points()?.iter().zip(strains) {
+            let force = young * area * e * p.weight * p.jacobian;
+            for (fi, ds) in f.iter_mut().zip(p.ds) {
+                for k in 0..3 {
+                    fi[k] += force * ds * p.tangent[k];
+                }
+            }
+        }
+        Ok(f)
+    }
+
+    /// Strain and stress at the integration points with an optional thermal strain: the
+    /// stress is E (ε − ε_th), the strain reported the total.
+    pub fn stresses_with(&self, young: f64, u: &[[f64; 3]], thermal: Option<&ThermalStrain>) -> Result<Vec<PointStress>> {
+        let strains = self.axial_strains(u)?;
+        let thermal = match thermal {
+            Some(t) => self.thermal_strains(t)?,
+            None => vec![0.0; strains.len()],
+        };
+        Ok(self
+            .points()?
+            .iter()
+            .zip(strains.into_iter().zip(thermal))
+            .map(|(p, (e, th))| PointStress { xyz: p.xyz, strain: uniaxial(e, p.tangent, 2.0), stress: uniaxial(young * (e - th), p.tangent, 1.0) })
             .collect())
     }
 
@@ -401,4 +440,26 @@ mod tests {
             assert!(err < 0.01 * span.powi(4), "{span}: {err}");
         }
     }
+
+    #[test]
+    fn thermal_strain_follows_the_axis() {
+        // A skew rod3 with orthotropic α: the free expansion along the member, u = ε_t (x − x₀)
+        // with ε_t = Σ αₖ tₖ² ΔT, is stress-free and balances the thermal load (K u = f_th).
+        let t = [2.0f64, 1.0, 2.0].map(|c| c / 3.0);
+        let xyz = [[0.0; 3], t.map(|c| 3.0 * c), t.map(|c| 1.5 * c)];
+        let truss = Truss::new(ElementKind::Rod3, &xyz).unwrap();
+        let alpha = [1.0e-5, 2.0e-5, 3.0e-5];
+        let th = ThermalStrain { alpha, delta: vec![40.0; 3] };
+        let eps: f64 = (0..3).map(|k| alpha[k] * t[k] * t[k]).sum::<f64>() * 40.0;
+        let u: Vec<[f64; 3]> = xyz.iter().map(|x| x.map(|c| eps * c)).collect();
+        for p in truss.stresses_with(2e5, &u, Some(&th)).unwrap() {
+            assert!(p.stress.iter().all(|s| s.abs() < 1e-9), "{:?}", p.stress);
+        }
+        let f = truss.thermal_load(2e5, 2.0, &th).unwrap();
+        let ku = truss.stiffness(2e5, 2.0).unwrap().mul_vec(&u.iter().flatten().copied().collect::<Vec<_>>());
+        for (a, b) in f.iter().flatten().zip(&ku) {
+            assert!((a - b).abs() < 1e-9 * 2e5 * 2.0 * eps, "{a} vs {b}");
+        }
+        assert!(truss.thermal_load(2e5, 2.0, &ThermalStrain { alpha, delta: vec![1.0] }).is_err());
+    }
 }