explicit: opt-in mass scaling to a target increment, the added mass reported

An explicit step may ask for mass scaling to a target time increment (on
a set, or the whole model): elements whose stable increment is below it
get the density that reaches it, with their lumped masses, frequency
bound and hourglass scales; adiabatic heating keeps the physical
density. The run reports the elements scaled, the added mass in all and
per section with its fraction, a warning, and the added mass's kinetic
energy in the energy account. Verified on a uniform bar scaled to twice
its increment (added mass 3×, increment doubled, a quarter of the
kinetic energy left to the physical mass). ROADMAP and parity note the
hyperelastic materials and mass scaling.

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

crates/explicit/src/elements.rs

@@ -547,6 +547,8 @@ fn hexahedron(sys: &System, h: &crate::system::Hex, m: &Motion, es: &mut Element
     let (d, _) = rotation::split(&g);
     let d = mechanical(d, h.thermal_rate * m.dt);
     let (rho, c) = current(mat, h.kinematics.volume, now.volume);
+    // Mass scaling raises the density the dynamics see (not the one adiabatic heating uses).
+    let (rho_s, c_s) = (rho * h.mass_factor, c / h.mass_factor.sqrt());
     let stiffer = match &mat.hyperelastic {
         Some(model) => {
             let mut f = h.kinematics.gradient(&e.u);
@@ -568,7 +570,7 @@ fn hexahedron(sys: &System, h: &crate::system::Hex, m: &Motion, es: &mut Element
     // geometry).
     let (gamma, scales) = if sys.nonlinear {
         let gamma = now.hourglass_vectors(&x_now);
-        let scales = hex_scales(&now, &gamma, &h.mass, mat, rho, hourglass);
+        let scales = hex_scales(&now, &gamma, &h.mass, mat, rho_s, hourglass);
         (gamma, scales)
     } else {
         (h.gamma, h.rest)
@@ -593,7 +595,7 @@ fn hexahedron(sys: &System, h: &crate::system::Hex, m: &Motion, es: &mut Element
             }
         }
     }
-    let element = Element { now: &now, mid: &mid, density: rho, wave_speed: c, length: 2.0 * c / (stiffer * scales.omega2).sqrt() };
+    let element = Element { now: &now, mid: &mid, density: rho_s, wave_speed: c_s, length: 2.0 * c_s / (stiffer * scales.omega2).sqrt() };
     let mut forces = element.finish(&es.point.stress, &e.v, m.bulk);
     forces.hourglass = Some(f);
     Ok(Outcome::Forces(forces))
@@ -632,6 +634,7 @@ fn tetrahedron(
     let thermal = t.thermal_rate * m.dt;
     let d = mechanical(d, thermal);
     let (rho, c) = current(mat, t.kinematics.volume, inc.now.volume);
+    let (rho_s, c_s) = (rho * t.mass_factor, c / t.mass_factor.sqrt());
     let stiffer = match &mat.hyperelastic {
         // The average nodal pressure of a decoupled material (Bonet & Burton 1998, at finite
         // strain): each node's pressure U′(J_a) from its tetrahedra's volume ratio, the element's
@@ -660,7 +663,7 @@ fn tetrahedron(
         return Ok(Outcome::Eroded(failure));
     }
     let omega2 = stiffer * if sys.nonlinear { frequency_bound(&inc.now, inc.now.gradient_norm(), mat.modulus, &t.mass) } else { t.omega2 };
-    let element = Element { now: &inc.now, mid: &inc.mid, density: rho, wave_speed: c, length: 2.0 * c / omega2.sqrt() };
+    let element = Element { now: &inc.now, mid: &inc.mid, density: rho_s, wave_speed: c_s, length: 2.0 * c_s / omega2.sqrt() };
     Ok(Outcome::Forces(element.finish(&es.point.stress, &inc.v, m.bulk)))
 }
 

crates/explicit/src/integrate.rs

@@ -201,10 +201,14 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res
     let mut critical: Vec<ControllingElement> = sys
         .hexes
         .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) })
+        .map(|h| (h.label, h.length, h.material, h.mass_factor))
+        .chain(sys.tets.iter().map(|t| (t.label, t.length, t.material, t.mass_factor)))
+        .chain(sys.rods.iter().map(|r| (r.label, r.length, r.material, 1.0)))
+        .map(|(label, length, m, f)| ControllingElement {
+            label,
+            // Mass scaling slows the element's waves by √f.
+            increment: length * f.sqrt() / 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);
@@ -411,6 +415,10 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res
         energy.external += external;
         energy.constraint += dt_n * constraint;
         energy.kinetic = 0.5 * kinetic;
+        if !sys.added_mass.is_empty() {
+            let per_node = v_now.as_chunks::<NODE_DOFS>().0.iter().zip(&sys.added_mass);
+            energy.added_kinetic = 0.5 * per_node.map(|(v, m)| m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2])).sum::<f64>();
+        }
         if !ws.contact_force.is_empty() {
             let work: f64 = du.iter().zip(ws.contact_force.iter().zip(&previous.contact_force)).map(|(d, (c, c0))| d * (c0 + c)).sum();
             energy.contact -= 0.5 * work;
@@ -483,6 +491,7 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res
         initial_kinetic,
         eroded,
         contacts: summaries,
+        mass_scaling: sys.mass_scaling.clone(),
         warnings,
     })
 }

crates/explicit/src/lib.rs

@@ -46,7 +46,7 @@ pub mod system;
 
 pub use elements::Failure;
 pub use integrate::solve_explicit;
-pub use result::{ContactSummary, ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History};
+pub use result::{ContactSummary, ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History, MassScalingSummary};
 
 #[derive(Debug, Clone, PartialEq, thiserror::Error)]
 pub enum Error {

crates/explicit/src/result.rs

@@ -25,6 +25,8 @@ pub struct Energies {
     /// Energy taken by contact: the work the model does against the contact forces (the penalty
     /// energy stored in the contact springs, plus what friction dissipates there).
     pub contact: f64,
+    /// The kinetic energy of the mass added by mass scaling (part of `kinetic`).
+    pub added_kinetic: f64,
     /// Friction dissipation, in contact and on rigid walls (part of `contact` and of
     /// `constraint`).
     pub friction: f64,
@@ -106,6 +108,19 @@ pub struct History {
     pub velocity: Vec<Vec<[f64; 3]>>,
 }
 
+/// What mass scaling added to a run.
+#[derive(Debug, Clone, Default, PartialEq, Serialize)]
+pub struct MassScalingSummary {
+    pub target_increment: f64,
+    /// Elements whose density was raised.
+    pub elements: usize,
+    pub added_mass: f64,
+    /// The added mass over the model's mass (with it).
+    pub fraction: f64,
+    /// The added mass per section, by name.
+    pub sections: Vec<(String, f64)>,
+}
+
 /// What a contact did over the run.
 #[derive(Debug, Clone, Default, PartialEq, Serialize)]
 pub struct ContactSummary {
@@ -158,5 +173,7 @@ pub struct ExplicitResult {
     pub eroded: Vec<ErodedElement>,
     /// Per contact of the model, in name order.
     pub contacts: Vec<ContactSummary>,
+    /// What mass scaling added, when the step asks for it.
+    pub mass_scaling: Option<MassScalingSummary>,
     pub warnings: Vec<String>,
 }

crates/explicit/src/system.rs

@@ -68,6 +68,8 @@ pub struct Hex {
     pub length: f64,
     /// The rate of its thermal strain, α ΔT / end time (every normal direction).
     pub thermal_rate: f64,
+    /// Its density's factor from mass scaling (1 without).
+    pub mass_factor: f64,
 }
 
 /// A nodal-pressure tetrahedron.
@@ -84,6 +86,8 @@ pub struct Tet {
     pub length: f64,
     /// The rate of its thermal strain, α ΔT / end time (every normal direction).
     pub thermal_rate: f64,
+    /// Its density's factor from mass scaling (1 without).
+    pub mass_factor: f64,
 }
 
 /// A fixed rigid plane.
@@ -120,6 +124,9 @@ pub struct System {
     /// 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)>,
+    /// What mass scaling added: per node, and per section (in name order).
+    pub added_mass: Vec<f64>,
+    pub mass_scaling: Option<crate::result::MassScalingSummary>,
     /// Prescribed velocity per degree of freedom (`None` = free).
     pub prescribed: Vec<Option<f64>>,
     pub walls: Vec<Wall>,
@@ -277,6 +284,8 @@ impl System {
             tets: Vec::new(),
             rods: Vec::new(),
             mass: vec![0.0; n],
+            added_mass: Vec::new(),
+            mass_scaling: None,
             discrete: Discrete::default(),
             discrete_increment: None,
             prescribed: vec![None; dofs],
@@ -322,7 +331,19 @@ impl System {
                     let rest = hex_scales(&k, &gamma, &mass, mat, rho, settings.hourglass);
                     sys.lump(&nodes, &mass);
                     let length = 2.0 * c / rest.omega2.sqrt();
-                    sys.hexes.push(Hex { label, index, nodes, material: mi, kinematics: k, gamma, mass, rest, length, thermal_rate: 0.0 });
+                    sys.hexes.push(Hex {
+                        label,
+                        index,
+                        nodes,
+                        material: mi,
+                        kinematics: k,
+                        gamma,
+                        mass,
+                        rest,
+                        length,
+                        thermal_rate: 0.0,
+                        mass_factor: 1.0,
+                    });
                 }
                 ElementKind::Tet4 => {
                     let nodes: [usize; 4] = idx.try_into().map_err(|_| invalid(format!("element {label} does not have 4 nodes")))?;
@@ -341,6 +362,7 @@ impl System {
                         omega2,
                         length: 2.0 * c / omega2.sqrt(),
                         thermal_rate: 0.0,
+                        mass_factor: 1.0,
                     });
                 }
                 ElementKind::Rod2 | ElementKind::Rod3 => {
@@ -369,6 +391,9 @@ impl System {
                 other => return Err(Error::Internal(format!("{} element kind", other.name()))),
             }
         }
+        if let Some(scaling) = &settings.mass_scaling {
+            sys.scale_mass(model, resolved, scaling, settings)?;
+        }
         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)?;
@@ -379,6 +404,90 @@ impl System {
         Ok(sys)
     }
 
+    /// Mass scaling to a target increment (see [`fieldcraft_model::MassScaling`]): each solid
+    /// element (of the set) whose stable increment times the step's scale is below the target
+    /// gets its density times (target / increment)², its lumped masses and frequency bound with
+    /// it; the added mass is kept per node and per section.
+    fn scale_mass(&mut self, model: &Model, resolved: &Resolved, scaling: &fieldcraft_model::MassScaling, settings: &ExplicitDynamic) -> Result<()> {
+        let chosen: Option<&[u32]> = match &scaling.set {
+            Some(s) => Some(resolved.sets.get(s).map(|r| r.elements.as_slice()).ok_or_else(|| invalid(format!("set `{s}` is not resolved")))?),
+            None => None,
+        };
+        let wanted = |label: u32| chosen.is_none_or(|l| l.binary_search(&label).is_ok());
+        let factor = |length: f64, material: usize, materials: &[Material]| -> f64 {
+            let c = materials.get(material).map_or(f64::INFINITY, |m| m.wave_speed);
+            let dt = settings.dt_scale * length / c;
+            if dt > 0.0 && dt < scaling.target_increment { (scaling.target_increment / dt).powi(2) } else { 1.0 }
+        };
+        let hourglass = fieldcraft_model::Hourglass { kind: settings.hourglass.kind, coefficient: settings.hourglass.coefficient };
+        let mut added = vec![0.0; self.labels.len()];
+        let mut by_element: Vec<(u32, f64)> = Vec::new();
+        for h in self.hexes.iter_mut().filter(|h| wanted(h.label)) {
+            let f = factor(h.length, h.material, &self.materials);
+            if f == 1.0 {
+                continue;
+            }
+            let extra = h.mass.map(|m| m * (f - 1.0));
+            for (n, x) in h.nodes.iter().zip(extra) {
+                if let Some(a) = added.get_mut(*n) {
+                    *a += x;
+                }
+            }
+            h.mass = h.mass.map(|m| m * f);
+            h.mass_factor = f;
+            if let Some(mat) = self.materials.get(h.material) {
+                h.rest = hex_scales(&h.kinematics, &h.gamma, &h.mass, mat, mat.density * f, hourglass);
+            }
+            by_element.push((h.label, extra.iter().sum()));
+        }
+        for t in self.tets.iter_mut().filter(|t| wanted(t.label)) {
+            let f = factor(t.length, t.material, &self.materials);
+            if f == 1.0 {
+                continue;
+            }
+            let extra = t.mass.map(|m| m * (f - 1.0));
+            for (n, x) in t.nodes.iter().zip(extra) {
+                if let Some(a) = added.get_mut(*n) {
+                    *a += x;
+                }
+            }
+            t.mass = t.mass.map(|m| m * f);
+            t.omega2 /= f;
+            t.mass_factor = f;
+            by_element.push((t.label, extra.iter().sum()));
+        }
+        for (m, a) in self.mass.iter_mut().zip(&added) {
+            *m += a;
+        }
+        let total_added: f64 = added.iter().sum();
+        let total: f64 = self.mass.iter().sum();
+        let mut sections: Vec<(String, f64)> = Vec::new();
+        for (name, section) in model.sections() {
+            let elements = resolved.sets.get(&section.set).map(|r| r.elements.as_slice()).unwrap_or_default();
+            let mass: f64 = by_element.iter().filter(|(l, _)| elements.binary_search(l).is_ok()).map(|(_, m)| m).sum();
+            if mass > 0.0 {
+                sections.push((name.clone(), mass));
+            }
+        }
+        if total_added > 0.0 {
+            self.warnings.push(format!(
+                "mass scaling to an increment of {:.3e} added {total_added:.4e} of mass ({:.2}% of the model's) to {} elements",
+                scaling.target_increment,
+                100.0 * total_added / total.max(f64::MIN_POSITIVE),
+                by_element.len()
+            ));
+        }
+        self.mass_scaling = Some(crate::result::MassScalingSummary {
+            target_increment: scaling.target_increment,
+            elements: by_element.len(),
+            added_mass: total_added,
+            fraction: total_added / total.max(f64::MIN_POSITIVE),
+            sections,
+        });
+        self.added_mass = added;
+        Ok(())
+    }
+
     /// Refuse what hyperelastic elements cannot take yet: small-displacement steps (the
     /// material is a finite-strain law) and thermal strain.
     fn check_hyperelastic(&self) -> Result<()> {

crates/explicit/tests/one_d.rs

@@ -305,3 +305,55 @@ fn whirling_mass_on_a_corotational_rod() {
     // A quarter turn: the mass ends near (0, L).
     assert!((l + end[0]).abs() < 1e-2 && (end[1] - l).abs() < 1e-2, "{end:?}");
 }
+
+/// Mass scaling to twice a uniform bar's stable increment quadruples its density: the added mass
+/// is three times the bar's, the stable increment doubles with the halved wave speed
+/// c = √(E/ρ), and the kinetic energy of the added mass is three quarters of the total.
+#[test]
+fn mass_scaling_slows_the_waves_as_the_added_density_says() {
+    use fieldcraft_explicit::{Limits, solve_explicit};
+    use fieldcraft_model::*;
+    let build = |scaled: Option<f64>| {
+        let mut m = Model::default();
+        let n = 40;
+        for i in 0..=n {
+            for (j, (y, z)) in [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)].iter().enumerate() {
+                m.set_node((4 * i + j + 1) as u32, [i as f64, *y, *z]).unwrap();
+            }
+        }
+        for i in 0..n {
+            let a = |j: usize| (4 * i + j + 1) as u32;
+            let b = |j: usize| (4 * (i + 1) + j + 1) as u32;
+            m.set_element(i as u32 + 1, ElementKind::Hex8, vec![a(0), b(0), b(1), a(1), a(3), b(3), b(2), a(2)]).unwrap();
+        }
+        m.add_set("all", Set { elements: (1..=n as u32).collect(), ..Set::default() }).unwrap();
+        m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: 1.0, poisson: 0.0 }), density: Some(1.0), ..Material::default() })
+            .unwrap();
+        let formulation = Formulation { integration: Integration::Reduced, ..Formulation::default() };
+        m.add_section("bar", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap();
+        let mut step = ExplicitDynamic::new(10.0);
+        step.initial_velocities = vec![InitialVelocity { set: "all".into(), velocity: [0.01, 0.0, 0.0] }];
+        step.mass_scaling = scaled.map(|t| MassScaling { target_increment: t, set: None });
+        m.add_step(Step::new("Run", Procedure::ExplicitDynamic(step)), None).unwrap();
+        solve_explicit(&m, None, &Limits::default()).unwrap()
+    };
+    let plain = build(None);
+    let scaled = build(Some(2.0 * plain.initial_increment));
+    let s = scaled.mass_scaling.as_ref().unwrap();
+    eprintln!(
+        "mass scaling: {s:?}; increments {} vs {}; Δt {:.4e} vs {:.4e}",
+        scaled.increments, plain.increments, scaled.initial_increment, plain.initial_increment
+    );
+    assert_eq!(s.elements, 40);
+    assert!((s.added_mass / plain.total_mass - 3.0).abs() < 1e-9 && (s.fraction - 0.75).abs() < 1e-9, "{s:?}");
+    assert_eq!(s.sections.len(), 1);
+    assert!((scaled.initial_increment / plain.initial_increment - 2.0).abs() < 1e-9);
+    assert!((scaled.total_mass / plain.total_mass - 4.0).abs() < 1e-9);
+    assert!((scaled.energy.added_kinetic / scaled.energy.kinetic - 0.75).abs() < 1e-9, "{:?}", scaled.energy);
+    assert!(scaled.warnings.iter().any(|w| w.contains("mass scaling")));
+    assert!(
+        scaled.controlling.first().is_some_and(|c| (0.9 * c.increment / scaled.initial_increment - 1.0).abs() < 1e-6),
+        "{:?}",
+        scaled.controlling
+    );
+}