forked from 1b-agent/fieldcraft-upstream
jyvikoandClaude Opus 5.5
a226b7dcac
At an inclined plane's nodes the velocity along its normal is removed every increment, from the components no prescribed velocity holds (planes meeting at a node orthonormalized, so one pass is exact); a prescribed motion off the plane is refused. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo
511 lines
24 KiB
Rust
511 lines
24 KiB
Rust
//! The time loop: central differences with variable increments, kinematic constraints, the
|
||
//! energy account and output.
|
||
//!
|
||
//! One increment (Belytschko, Liu, Moran, Elkhodary 2014, Box 6.1, with Δt^n = ½(Δt^{n−½} +
|
||
//! Δt^{n+½})):
|
||
//! 1. u^n = u^{n−1} + Δt^{n−½} v^{n−½}; t^n = t^{n−1} + Δt^{n−½};
|
||
//! 2. internal forces f^n from u^n (stresses, hourglass stiffness) and v^{n−½} (hourglass and
|
||
//! bulk viscosity), and the next increment Δt^{n+½} = scale · min_e Δt_e;
|
||
//! 3. the contact forces c^n at u^n ([`crate::contact`], with the stable increment Δt^{n+½});
|
||
//! v^{n+½} = v^{n−½} + Δt^n M⁻¹(f_ext + c^n − f^n), then the constraints: rigid walls set the
|
||
//! normal velocity of a node that would cross its wall so it lands on it (a kinematic contact
|
||
//! condition: T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier, W. Kanoknukulchai,
|
||
//! *A finite element method for a class of contact-impact problems*, CMAME 8(3):249–276,
|
||
//! 1976), with Coulomb friction as a tangential impulse at most μ times the normal one (the
|
||
//! tangential velocity stops, or loses μ Δv_n when it slides), prescribed degrees of
|
||
//! freedom take their velocity, and a symmetry plane not normal to an axis removes its nodes'
|
||
//! velocity along its normal (from their free components: the lumped mass is the same on a
|
||
//! node's three translations, so this is the mass-weighted projection onto the constraint). The change a constraint makes is the impulse of its reaction,
|
||
//! r^n = M Δv / Δt^n;
|
||
//! 4. energies: the works of the internal, contact and external forces by the trapezoidal rule,
|
||
//! W += ½ Δu · (f^{n−1} + f^n) (exactly the strain energy for a linear elastic material),
|
||
//! the constraints' work Δt^n r^n · ½(v^{n−½} + v^{n+½}) (exactly the kinetic energy their
|
||
//! impulses add or remove), and the kinetic energy of v^n = v^{n−½} + ½Δt^{n−½} a^n.
|
||
//!
|
||
//! The run starts from u = 0 and the initial velocities with the half increment
|
||
//! v^{½} = v^0 + ½Δt^{½} a^0.
|
||
|
||
/// Increments between progress reports (each takes a lock).
|
||
const PROGRESS_EVERY: u64 = 64;
|
||
|
||
use fieldcraft_model::{ExplicitDynamic, Model, OutputVariable, Procedure, Step};
|
||
|
||
use crate::contact::{Evaluation, Interface};
|
||
use crate::elements::{Motion, State, Workspace, internal_forces};
|
||
use crate::result::{ContactSummary, ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History};
|
||
use crate::system::{NODE_DOFS, System};
|
||
use crate::{Error, Limits, Result};
|
||
|
||
/// How many elements the result names as controlling the time increment.
|
||
pub const CONTROLLING: usize = 10;
|
||
|
||
/// Relative energy balance error above which a run warns.
|
||
pub const BALANCE_TOLERANCE: f64 = 0.01;
|
||
|
||
/// Hourglass energy as a fraction of the largest internal energy of the run above which it warns
|
||
/// (the internal energy may fall back to zero, as in a rebound).
|
||
pub const HOURGLASS_TOLERANCE: f64 = 0.1;
|
||
|
||
/// The step named `name`, or the model's only explicit dynamic step.
|
||
fn select_step<'m>(model: &'m Model, name: Option<&str>) -> Result<(&'m Step, &'m ExplicitDynamic)> {
|
||
let step = model.select_step(name, Procedure::is_dynamic, "explicit dynamic")?;
|
||
match &step.procedure {
|
||
Procedure::ExplicitDynamic(e) => Ok((step, e)),
|
||
_ => Err(Error::Internal("an explicit dynamic step".into())),
|
||
}
|
||
}
|
||
|
||
/// Times at a fixed interval: whether `t` has reached the next one.
|
||
struct Schedule {
|
||
interval: f64,
|
||
next: f64,
|
||
}
|
||
|
||
impl Schedule {
|
||
fn due(&mut self, t: f64) -> bool {
|
||
if t < self.next * (1.0 - 1e-12) {
|
||
return false;
|
||
}
|
||
self.next = ((t / self.interval).floor() + 1.0) * self.interval;
|
||
true
|
||
}
|
||
}
|
||
|
||
/// What the run writes.
|
||
struct Output {
|
||
displacement: bool,
|
||
stress: bool,
|
||
strain: bool,
|
||
}
|
||
|
||
/// Velocity constraints over an increment of length `dt` (rigid walls when dt > 0, then
|
||
/// prescribed degrees of freedom, then inclined symmetry planes); returns the kinetic energy the walls' friction took (`mass`
|
||
/// per degree of freedom).
|
||
fn constrain(sys: &System, u: &[f64], v: &mut [f64], dt: f64, mass: &[f64]) -> f64 {
|
||
let mut dissipated = 0.0;
|
||
if dt > 0.0 {
|
||
for w in &sys.walls {
|
||
for &node in &w.nodes {
|
||
let x = sys.xyz.get(node).copied().unwrap_or_default();
|
||
let at = |k: usize| NODE_DOFS * node + k;
|
||
let gap: f64 = (0..3).map(|k| (x[k] + u.get(at(k)).copied().unwrap_or(0.0) - w.point[k]) * w.normal[k]).sum();
|
||
let normal: f64 = (0..3).map(|k| v.get(at(k)).copied().unwrap_or(0.0) * w.normal[k]).sum();
|
||
if gap + dt * normal < 0.0 {
|
||
let change = -gap / dt - normal;
|
||
for k in 0..3 {
|
||
if let Some(c) = v.get_mut(at(k)) {
|
||
*c += change * w.normal[k];
|
||
}
|
||
}
|
||
if !w.friction.is_frictionless() {
|
||
dissipated += wall_friction(v, node, w, change, mass.get(at(0)).copied().unwrap_or(0.0));
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
for (c, p) in v.iter_mut().zip(&sys.prescribed) {
|
||
if let Some(p) = p {
|
||
*c = *p;
|
||
}
|
||
}
|
||
for row in &sys.planes {
|
||
row.project(v);
|
||
}
|
||
dissipated
|
||
}
|
||
|
||
/// Coulomb friction of a node the wall `w` stopped with a normal velocity change `change`: its
|
||
/// tangential velocity changes by at most μ · change (stick when that stops it); returns the
|
||
/// kinetic energy taken, ½ m (|v_t|² − |v_t′|²).
|
||
fn wall_friction(v: &mut [f64], node: usize, w: &crate::system::Wall, change: f64, mass: f64) -> f64 {
|
||
let at = |k: usize| NODE_DOFS * node + k;
|
||
let vel = [0, 1, 2].map(|k| v.get(at(k)).copied().unwrap_or(0.0));
|
||
let vn: f64 = (0..3).map(|k| vel[k] * w.normal[k]).sum();
|
||
let tangential = [0, 1, 2].map(|k| vel[k] - vn * w.normal[k]);
|
||
let slide = tangential.iter().map(|c| c * c).sum::<f64>().sqrt();
|
||
if slide == 0.0 {
|
||
return 0.0;
|
||
}
|
||
let cap = w.friction.coefficient(slide) * change;
|
||
let keep = if slide <= cap { 0.0 } else { 1.0 - cap / slide };
|
||
for k in 0..3 {
|
||
if let Some(c) = v.get_mut(at(k)) {
|
||
*c = vn * w.normal[k] + keep * tangential[k];
|
||
}
|
||
}
|
||
0.5 * mass * slide * slide * (1.0 - keep * keep)
|
||
}
|
||
|
||
/// v_out = v_in + dt M⁻¹ (f_ext + f_contact − f_int).
|
||
fn accelerate(sys: &System, ws: &Workspace, mass: &[f64], v_in: &[f64], dt: f64, v_out: &mut [f64]) {
|
||
let internal = ws.stress_force.iter().zip(&ws.hourglass_force).zip(&ws.viscous_force);
|
||
for ((((out, v), m), e), ((s, h), q)) in v_out.iter_mut().zip(v_in).zip(mass).zip(&sys.external).zip(internal) {
|
||
*out = v + dt * ((e - (s + h + q)) / m);
|
||
}
|
||
// Contact adds its forces (a separate pass, so a model without contact runs the same sums).
|
||
if !ws.contact_force.is_empty() {
|
||
for ((out, m), c) in v_out.iter_mut().zip(mass).zip(&ws.contact_force) {
|
||
*out += dt * (c / m);
|
||
}
|
||
}
|
||
}
|
||
|
||
/// The contacts' forces at `u` (with the increment `du` over `dt_prev` that led there, and the
|
||
/// stable increment `dt`) into `ws.contact_force`, and what each found.
|
||
#[allow(clippy::too_many_arguments)]
|
||
fn contact_forces(
|
||
contacts: &mut [Interface],
|
||
sys: &System,
|
||
u: &[f64],
|
||
du: &[f64],
|
||
dt_prev: f64,
|
||
dt: f64,
|
||
mass: &[f64],
|
||
state: &State,
|
||
ws: &mut Workspace,
|
||
) -> Vec<Evaluation> {
|
||
ws.contact_force.fill(0.0);
|
||
contacts.iter_mut().map(|c| c.forces(&sys.xyz, u, du, dt_prev, dt, mass, state, &mut ws.contact_force)).collect()
|
||
}
|
||
|
||
/// A contact's record over the run.
|
||
fn track(summary: &mut ContactSummary, e: &Evaluation, t: f64) {
|
||
if e.active > 0 {
|
||
summary.first_contact.get_or_insert(t);
|
||
summary.last_contact = Some(t);
|
||
}
|
||
if e.normal_force > summary.peak_force {
|
||
summary.peak_force = e.normal_force;
|
||
summary.peak_time = t;
|
||
}
|
||
summary.max_penetration = summary.max_penetration.max(e.penetration);
|
||
summary.friction += e.dissipation;
|
||
}
|
||
|
||
/// The constraint reactions M (v_constrained − v_free) / Δt, into `r`.
|
||
fn impulses(mass: &[f64], free: &[f64], constrained: &[f64], dt: f64, r: &mut [f64]) {
|
||
for (r, (m, (f, c))) in r.iter_mut().zip(mass.iter().zip(free.iter().zip(constrained))) {
|
||
*r = if c == f { 0.0 } else { m * (c - f) / dt };
|
||
}
|
||
}
|
||
|
||
fn per_node(x: &[f64]) -> Vec<[f64; 3]> {
|
||
x.as_chunks::<NODE_DOFS>().0.to_vec()
|
||
}
|
||
|
||
/// Solve an explicit dynamic step: the step named `step`, or the model's only one.
|
||
pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Result<ExplicitResult> {
|
||
let resolved = model.resolve()?;
|
||
let (step, settings) = select_step(model, step)?;
|
||
let mut sys = System::new(model, &resolved, step, settings, limits)?;
|
||
let end = settings.end_time;
|
||
let scale = settings.dt_scale;
|
||
let bulk = settings.bulk_viscosity;
|
||
let max_increments = settings.max_increments.min(limits.max_increments);
|
||
let mut warnings = sys.warnings.clone();
|
||
|
||
// The controlling elements: the smallest critical increments at rest.
|
||
let mut critical: Vec<ControllingElement> = sys
|
||
.hexes
|
||
.iter()
|
||
.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);
|
||
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);
|
||
if estimate > max_increments as f64 {
|
||
let first = critical.first().map_or(0, |c| c.label);
|
||
return Err(Error::Limit(format!(
|
||
"the step needs about {estimate:.3e} increments of {:.3e} (element {first} controls it), more than the cap of {max_increments}",
|
||
scale * smallest
|
||
)));
|
||
}
|
||
|
||
// Output sizes.
|
||
let output = Output {
|
||
displacement: step.output.contains(&OutputVariable::Displacement),
|
||
stress: step.output.contains(&OutputVariable::Stress),
|
||
strain: step.output.contains(&OutputVariable::Strain),
|
||
};
|
||
let (nodes, elements) = (sys.labels.len(), sys.elements.len());
|
||
let frames = (end / settings.output_interval()).ceil() + 2.0;
|
||
let samples = (end / settings.history_interval()).ceil() + 2.0;
|
||
let per_frame = (if output.displacement { 6 * nodes } else { 0 }) + 6 * elements * (usize::from(output.stress) + usize::from(output.strain));
|
||
let per_sample = 16 + 6 * sys.history_nodes.len();
|
||
let values = frames * per_frame as f64 + samples * per_sample as f64;
|
||
if values > limits.max_output_values as f64 {
|
||
return Err(Error::Limit(format!(
|
||
"the output would hold {values:.3e} numbers; the limit is {} (lengthen the output or history interval, or request less output)",
|
||
limits.max_output_values
|
||
)));
|
||
}
|
||
|
||
// Point masses count with their mean over the three directions.
|
||
let total_mass: f64 = sys.mass.iter().sum::<f64>() + sys.discrete.mass.iter().sum::<f64>() / 3.0;
|
||
let mut ws = Workspace::new(&sys);
|
||
let mut previous = Workspace::new(&sys);
|
||
let mut contacts = std::mem::take(&mut sys.contacts);
|
||
let mut summaries: Vec<ContactSummary> =
|
||
contacts.iter().map(|c| ContactSummary { name: c.name.clone(), adjusted: c.adjusted.clone(), ..ContactSummary::default() }).collect();
|
||
let mut state = State::new(&sys);
|
||
let mut eroded: Vec<ErodedElement> = Vec::new();
|
||
let mut u = vec![0.0; n];
|
||
let mut du = vec![0.0; n];
|
||
let v0 = sys.initial_velocity.clone();
|
||
let initial_kinetic = 0.5 * v0.iter().zip(&mass).map(|(v, m)| m * v * v).sum::<f64>();
|
||
let inelastic = sys.materials.iter().any(|m| m.j2.as_ref().is_some_and(|j| j.flow.is_some()));
|
||
|
||
// Start: forces at rest with the initial velocities, then the half increment.
|
||
let start = Motion { u: &u, du: &du, v: &v0, dt: 0.0, time: 0.0, bulk: &bulk };
|
||
let crit = internal_forces(&sys, &start, &mut state, &mut ws)?;
|
||
let mut dt_half = (scale * crit).min(end);
|
||
let initial_increment = dt_half;
|
||
let mut found = contact_forces(&mut contacts, &sys, &u, &du, 0.0, scale * crit, &mass, &state, &mut ws);
|
||
for (summary, e) in summaries.iter_mut().zip(&found) {
|
||
track(summary, e, 0.0);
|
||
}
|
||
let mut v_half = vec![0.0; n];
|
||
accelerate(&sys, &ws, &mass, &v0, 0.5 * dt_half, &mut v_half);
|
||
let mut free = v_half.clone();
|
||
let wall_friction = constrain(&sys, &u, &mut v_half, dt_half, &mass);
|
||
let mut reaction = vec![0.0; n];
|
||
impulses(&mass, &free, &v_half, 0.5 * dt_half, &mut reaction);
|
||
// The constraints' work is that of their impulses on the mean of the velocities before
|
||
// and after the increment, Δt r · ½(v^{n−½} + v^{n+½}): exact for the kinetic energy a
|
||
// kinematic constraint removes or adds (a node stopped by a wall loses ½ m v²).
|
||
let impulse_work = |r: &[f64], before: &[f64], after: &[f64], dt: f64| -> f64 {
|
||
dt * r.iter().zip(before.iter().zip(after)).map(|(r, (a, b))| r * 0.5 * (a + b)).sum::<f64>()
|
||
};
|
||
let mut energy = Energies {
|
||
kinetic: initial_kinetic,
|
||
constraint: impulse_work(&reaction, &v0, &v_half, 0.5 * dt_half),
|
||
friction: wall_friction,
|
||
..Energies::default()
|
||
};
|
||
|
||
let mut history = History { nodes: sys.history_nodes.iter().filter_map(|i| sys.labels.get(*i).copied()).collect(), ..History::default() };
|
||
let mut out_frames = Vec::new();
|
||
let sample = |history: &mut History, t: f64, dt: f64, e: Energies, u: &[f64], v: &[f64], r: &[f64], found: &[Evaluation]| {
|
||
history.time.push(t);
|
||
history.contact_force.push(found.iter().map(|e| e.normal_force).collect());
|
||
history.increment.push(dt);
|
||
history.energies.push(e);
|
||
let sum = |x: &[f64], weight: &dyn Fn(usize, usize) -> f64| {
|
||
let mut s = [0.0; 3];
|
||
for (node, c) in x.as_chunks::<NODE_DOFS>().0.iter().enumerate() {
|
||
for k in 0..3 {
|
||
s[k] += weight(node, k) * c[k];
|
||
}
|
||
}
|
||
s
|
||
};
|
||
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));
|
||
};
|
||
let frame = |t: f64, u: &[f64], v: &[f64], state: &State| {
|
||
let each = |on: bool, f: &dyn Fn(&crate::elements::ElementState) -> f64| -> Vec<f64> {
|
||
if on { state.elements.iter().map(f).collect() } else { Vec::new() }
|
||
};
|
||
let inelastic = output.stress && inelastic;
|
||
Frame {
|
||
time: t,
|
||
displacement: if output.displacement { per_node(u) } else { Vec::new() },
|
||
velocity: if output.displacement { per_node(v) } else { Vec::new() },
|
||
stress: if output.stress { state.elements.iter().map(|e| e.point.stress).collect() } else { Vec::new() },
|
||
strain: if output.strain { state.elements.iter().map(|e| e.strain).collect() } else { Vec::new() },
|
||
plastic_strain: each(inelastic, &|e| e.point.plastic_strain),
|
||
damage: each(inelastic, &|e| e.point.damage),
|
||
temperature: each(inelastic, &|e| e.point.temperature),
|
||
failed: if output.stress { state.elements.iter().map(|e| e.failed.is_some()).collect() } else { Vec::new() },
|
||
}
|
||
};
|
||
let plastic = |state: &State| state.elements.iter().map(|e| e.plastic_work).sum::<f64>();
|
||
let mut histories = Schedule { interval: settings.history_interval(), next: settings.history_interval() };
|
||
let mut frames_due = Schedule { interval: settings.output_interval(), next: settings.output_interval() };
|
||
sample(&mut history, 0.0, dt_half, energy, &u, &v0, &reaction, &found);
|
||
out_frames.push(frame(0.0, &u, &v0, &state));
|
||
|
||
let mut t = 0.0;
|
||
let mut increments = 0u64;
|
||
let mut peak_internal: f64 = 0.0;
|
||
let mut v_now = v0;
|
||
let mut completed = true;
|
||
let mut v_new = vec![0.0; n];
|
||
fieldcraft_model::progress::stage("Time", Some(end), "");
|
||
loop {
|
||
if increments.is_multiple_of(PROGRESS_EVERY) {
|
||
fieldcraft_model::progress::count(t);
|
||
fieldcraft_model::progress::check()?;
|
||
}
|
||
if increments >= max_increments {
|
||
completed = false;
|
||
warnings.push(format!("stopped at time {t:.6e} after {increments} increments (the cap); the step ends at {end:.6e}"));
|
||
break;
|
||
}
|
||
for ((x, dx), v) in u.iter_mut().zip(du.iter_mut()).zip(&v_half) {
|
||
*dx = dt_half * v;
|
||
*x += *dx;
|
||
}
|
||
t += dt_half;
|
||
increments += 1;
|
||
std::mem::swap(&mut ws, &mut previous);
|
||
let motion = Motion { u: &u, du: &du, v: &v_half, dt: dt_half, time: t, bulk: &bulk };
|
||
let crit = internal_forces(&sys, &motion, &mut state, &mut ws)?;
|
||
for &(index, failure) in &state.eroded {
|
||
let element = state.elements.get(index);
|
||
let internal = element.map_or(0.0, |e| e.energy);
|
||
energy.eroded += internal;
|
||
eroded.push(ErodedElement { label: sys.elements.get(index).copied().unwrap_or(0), time: t, failure, internal_energy: internal });
|
||
}
|
||
let remaining = end - t;
|
||
let last = remaining <= 1e-12 * end;
|
||
let dt_next = if last { 0.0 } else { (scale * crit).min(remaining) };
|
||
let dt_n = 0.5 * (dt_half + dt_next);
|
||
if !contacts.is_empty() {
|
||
found = contact_forces(&mut contacts, &sys, &u, &du, dt_half, scale * crit, &mass, &state, &mut ws);
|
||
for (summary, e) in summaries.iter_mut().zip(&found) {
|
||
track(summary, e, t);
|
||
energy.friction += e.dissipation;
|
||
}
|
||
}
|
||
accelerate(&sys, &ws, &mass, &v_half, dt_n, &mut v_new);
|
||
free.clone_from(&v_new);
|
||
energy.friction += constrain(&sys, &u, &mut v_new, dt_next, &mass);
|
||
impulses(&mass, &free, &v_new, dt_n, &mut reaction);
|
||
// One pass over the degrees of freedom: v^n = v^{n−½} + ½Δt^{n−½} a^n with the constrained
|
||
// a^n = (v^{n+½} − v^{n−½}) / Δt^n, and the works (trapezoidal over the displacement for
|
||
// the internal and external forces).
|
||
let weight = 0.5 * dt_half / dt_n;
|
||
let [mut internal, mut hourglass, mut viscous, mut external, mut constraint, mut kinetic] = [0.0; 6];
|
||
let forces = ws.stress_force.iter().zip(&ws.hourglass_force).zip(&ws.viscous_force);
|
||
let before = previous.stress_force.iter().zip(&previous.hourglass_force).zip(&previous.viscous_force);
|
||
let per_dof = v_now.iter_mut().zip(&v_half).zip(&v_new).zip(&mass).zip(&du).zip(&sys.external).zip(&reaction);
|
||
for ((((((((now, vh), vn), m), dx), e), r), ((s, h), q)), ((s0, h0), q0)) in per_dof.zip(forces).zip(before) {
|
||
*now = vh + weight * (vn - vh);
|
||
internal += dx * (s0 + s);
|
||
hourglass += dx * (h0 + h);
|
||
viscous += dx * (q0 + q);
|
||
external += dx * e;
|
||
constraint += r * 0.5 * (vh + vn);
|
||
kinetic += m * *now * *now;
|
||
}
|
||
energy.internal += 0.5 * internal;
|
||
energy.hourglass += 0.5 * hourglass;
|
||
energy.viscous += 0.5 * viscous;
|
||
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;
|
||
}
|
||
peak_internal = peak_internal.max(energy.internal.abs());
|
||
energy.balance = energy.kinetic + energy.internal + energy.hourglass + energy.viscous + energy.contact
|
||
- energy.external
|
||
- energy.constraint
|
||
- initial_kinetic;
|
||
// Every element eroded leaves an infinite critical increment: the nodes fly freely.
|
||
if !(energy.balance.is_finite() && !crit.is_nan()) {
|
||
return Err(Error::Unstable(format!("the solution is not finite at time {t:.6e} (increment {increments})")));
|
||
}
|
||
if histories.due(t) || last {
|
||
energy.plastic = plastic(&state);
|
||
sample(&mut history, t, dt_half, energy, &u, &v_now, &reaction, &found);
|
||
}
|
||
if frames_due.due(t) || last {
|
||
out_frames.push(frame(t, &u, &v_now, &state));
|
||
}
|
||
std::mem::swap(&mut v_half, &mut v_new);
|
||
dt_half = dt_next;
|
||
if last {
|
||
break;
|
||
}
|
||
}
|
||
if !completed {
|
||
if history.time.last() != Some(&t) {
|
||
sample(&mut history, t, dt_half, energy, &u, &v_now, &reaction, &found);
|
||
}
|
||
if out_frames.last().map(|f| f.time) != Some(t) {
|
||
out_frames.push(frame(t, &u, &v_now, &state));
|
||
}
|
||
}
|
||
energy.plastic = plastic(&state);
|
||
if !eroded.is_empty() {
|
||
warnings.push(format!(
|
||
"{} elements eroded (their material failed); their internal energy at erosion was {:.4e}",
|
||
eroded.len(),
|
||
energy.eroded
|
||
));
|
||
}
|
||
|
||
let error = energy.relative_error(initial_kinetic);
|
||
if error > BALANCE_TOLERANCE {
|
||
warnings.push(format!("the energy balance is off by {:.2}% of the largest energy (kinetic, internal or work)", 100.0 * error));
|
||
}
|
||
// Round-off hourglass energy in a model at rest or in rigid motion is not worth a warning.
|
||
let floor =
|
||
1e-9 * [initial_kinetic, energy.kinetic, peak_internal, energy.external.abs(), energy.constraint.abs()].into_iter().fold(0.0, f64::max);
|
||
if energy.hourglass.abs() > HOURGLASS_TOLERANCE * peak_internal && energy.hourglass.abs() > floor {
|
||
warnings.push(format!(
|
||
"hourglass energy is {:.1}% of the largest internal energy; refine the mesh or change the hourglass control",
|
||
100.0 * energy.hourglass.abs() / peak_internal.max(f64::MIN_POSITIVE)
|
||
));
|
||
}
|
||
Ok(ExplicitResult {
|
||
step: step.name.clone(),
|
||
nodes: sys.labels.clone(),
|
||
elements: sys.elements.clone(),
|
||
frames: out_frames,
|
||
history,
|
||
controlling: critical,
|
||
initial_increment,
|
||
increments,
|
||
time: t,
|
||
completed,
|
||
total_mass,
|
||
energy,
|
||
initial_kinetic,
|
||
eroded,
|
||
contacts: summaries,
|
||
mass_scaling: sys.mass_scaling.clone(),
|
||
warnings,
|
||
})
|
||
}
|