forked from 1b-agent/fieldcraft-upstream
jyvikoandClaude Opus 5.5
bc9fe2191a
A nonlinear_static procedure (nonlinear geometry, initial/min/max increments, increment and iteration caps, force/displacement/energy tolerances, line search), validated with caps against hostile values. Materials gain an explicit hyperelastic option (decoupled compressible neo-Hookean on the isotropic constants). Abaqus writes *STEP, NLGEOM, INC / *STATIC with the increments and *HYPERELASTIC, NEO HOOKE; LS-DYNA writes the ramped implicit run with full Newton and *CONTROL_IMPLICIT_AUTO; Nastran refuses SOL 400 by name, as LS-DYNA does hyperelastic materials in nonlinear steps. Explicit dynamics refuses hyperelastic materials by name; solve.run names the solver still to come. Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm
599 lines
29 KiB
Rust
599 lines
29 KiB
Rust
//! From the resolved model to the explicit system: elements (solids and trusses) with their
|
||
//! materials and stable lengths, lumped masses with the point masses, springs and dashpots
|
||
//! ([`crate::discrete`]), prescribed velocities, rigid walls, initial velocities and the load
|
||
//! vector of an explicit dynamic step.
|
||
|
||
use fieldcraft_fem::material::Elasticity;
|
||
use fieldcraft_fem::plasticity::J2;
|
||
use fieldcraft_fem::uniform::{Hex8, Tet4, modulus_bound};
|
||
use fieldcraft_fem::{DofKind, DofLayout, DofSet, face};
|
||
use fieldcraft_model::{
|
||
Elastic, ElementKind, ExplicitDynamic, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step, StressRate, Target,
|
||
};
|
||
|
||
use crate::discrete::{Discrete, Rod};
|
||
use crate::elements::{HexScales, add_at, frequency_bound, hex_scales};
|
||
use crate::{Error, Result};
|
||
|
||
/// Degrees of freedom per node: the three translations, node-major (see [`System`]).
|
||
pub const NODE_DOFS: usize = 3;
|
||
|
||
/// Most nodes whose histories a step records.
|
||
pub const MAX_HISTORY_NODES: usize = 1000;
|
||
|
||
/// A material as the elements use it.
|
||
#[derive(Debug, Clone)]
|
||
pub struct Material {
|
||
pub name: String,
|
||
pub elasticity: Elasticity,
|
||
pub density: f64,
|
||
/// Dilatational wave speed √(M/ρ), with M the largest normal stiffness (λ + 2μ when
|
||
/// isotropic).
|
||
pub wave_speed: f64,
|
||
/// The bound of [`fieldcraft_fem::uniform::modulus_bound`].
|
||
pub modulus: f64,
|
||
/// The largest normal stiffness M (λ + 2μ when isotropic): c = √(M/ρ).
|
||
pub stiffest: f64,
|
||
/// The hypoelastic(–plastic) update of an isotropic material; `None` for orthotropic
|
||
/// elasticity, updated with [`Material::elasticity`].
|
||
pub j2: Option<J2>,
|
||
/// The equivalent plastic strain at which elements erode.
|
||
pub erosion: Option<f64>,
|
||
/// The (isotropic) expansion coefficient and reference temperature, when the material
|
||
/// expands.
|
||
pub expansion: Option<(f64, f64)>,
|
||
}
|
||
|
||
/// A one-point hexahedron.
|
||
#[derive(Debug, Clone)]
|
||
pub struct Hex {
|
||
pub label: u32,
|
||
/// Position in [`System::elements`].
|
||
pub index: usize,
|
||
pub nodes: [usize; 8],
|
||
pub material: usize,
|
||
pub kinematics: Hex8,
|
||
pub gamma: [[f64; 8]; 4],
|
||
/// Row-sum lumped masses of the element's nodes.
|
||
pub mass: [f64; 8],
|
||
/// Size, hourglass coefficient and frequency bound of the undeformed element.
|
||
pub rest: HexScales,
|
||
/// The element's stable length: its critical increment without bulk viscosity is
|
||
/// `length / wave_speed`.
|
||
pub length: f64,
|
||
/// The rate of its thermal strain, α ΔT / end time (every normal direction).
|
||
pub thermal_rate: f64,
|
||
}
|
||
|
||
/// A nodal-pressure tetrahedron.
|
||
#[derive(Debug, Clone)]
|
||
pub struct Tet {
|
||
pub label: u32,
|
||
pub index: usize,
|
||
pub nodes: [usize; 4],
|
||
pub material: usize,
|
||
pub kinematics: Tet4,
|
||
pub mass: [f64; 4],
|
||
/// The frequency bound ω² of the undeformed element.
|
||
pub omega2: f64,
|
||
pub length: f64,
|
||
/// The rate of its thermal strain, α ΔT / end time (every normal direction).
|
||
pub thermal_rate: f64,
|
||
}
|
||
|
||
/// A fixed, frictionless rigid plane.
|
||
#[derive(Debug, Clone)]
|
||
pub struct Wall {
|
||
pub point: [f64; 3],
|
||
/// Unit normal, towards the side the nodes stay on.
|
||
pub normal: [f64; 3],
|
||
pub nodes: Vec<usize>,
|
||
}
|
||
|
||
/// The discrete system of an explicit dynamic step.
|
||
///
|
||
/// Its elements are solids, so every node carries the three translations and the layout is
|
||
/// node-major: degree of freedom 3 n + axis, which the integrator's per-node loops rely on
|
||
/// (checked when the system is built).
|
||
#[derive(Debug, Clone)]
|
||
pub struct System {
|
||
/// Node labels with degrees of freedom (those of the elements), sorted.
|
||
pub labels: Vec<u32>,
|
||
pub layout: DofLayout,
|
||
pub xyz: Vec<[f64; 3]>,
|
||
/// Element labels, sorted.
|
||
pub elements: Vec<u32>,
|
||
pub materials: Vec<Material>,
|
||
pub hexes: Vec<Hex>,
|
||
pub tets: Vec<Tet>,
|
||
pub rods: Vec<Rod>,
|
||
/// Row-sum lumped mass per node.
|
||
pub mass: Vec<f64>,
|
||
/// Point masses, springs and dashpots.
|
||
pub discrete: Discrete,
|
||
/// The springs' and dashpots' stable increment on the lumped mass, and a link that sets it,
|
||
/// when there are any.
|
||
pub discrete_increment: Option<(f64, String)>,
|
||
/// Prescribed velocity per degree of freedom (`None` = free).
|
||
pub prescribed: Vec<Option<f64>>,
|
||
pub walls: Vec<Wall>,
|
||
/// Velocity per degree of freedom at the start (prescribed velocities included).
|
||
pub initial_velocity: Vec<f64>,
|
||
/// Constant external force per degree of freedom.
|
||
pub external: Vec<f64>,
|
||
/// Nodes whose histories are recorded.
|
||
pub history_nodes: Vec<usize>,
|
||
pub hourglass: HourglassKind,
|
||
pub hourglass_coefficient: f64,
|
||
/// Large deformation (updated Lagrangian) or the undeformed mesh's gradients.
|
||
pub nonlinear: bool,
|
||
pub stress_rate: StressRate,
|
||
pub warnings: Vec<String>,
|
||
}
|
||
|
||
fn invalid(m: impl Into<String>) -> Error {
|
||
Error::Invalid(m.into())
|
||
}
|
||
|
||
/// A material's constants for the explicit elements.
|
||
fn material(model: &Model, name: &str) -> Result<Material> {
|
||
let m = model.material(name).ok_or_else(|| invalid(format!("material `{name}` is missing")))?;
|
||
let elastic = m.elastic.as_ref().ok_or_else(|| invalid(format!("material `{name}` has no elastic constants")))?;
|
||
if let Some(h) = m.hyperelastic {
|
||
return Err(Error::NotSupported(format!("{} hyperelastic material `{name}` in explicit dynamics", h.name())));
|
||
}
|
||
let density = m
|
||
.density
|
||
.filter(|d| *d > 0.0 && d.is_finite())
|
||
.ok_or_else(|| invalid(format!("material `{name}` has no density, and explicit dynamics needs one")))?;
|
||
let d = fieldcraft_fem::elasticity(elastic).map_err(|e| invalid(format!("material `{name}`: {e}")))?;
|
||
let stiffest = d[0][0].max(d[1][1]).max(d[2][2]);
|
||
// The model keeps plasticity, damage and erosion to isotropic materials.
|
||
let j2 = match *elastic {
|
||
Elastic::Isotropic { .. } => {
|
||
Some(J2::from_model(elastic, m.plastic.as_ref(), m.damage.as_ref()).map_err(|e| invalid(format!("material `{name}`: {e}")))?)
|
||
}
|
||
Elastic::EngineeringConstants { .. } => None,
|
||
};
|
||
let expansion = match m.thermal.and_then(|t| t.expansion.map(|a| (a.axes(), t.reference_temperature))) {
|
||
None => None,
|
||
Some(([a, b, c], t_ref)) if a == b && b == c => Some((a, t_ref)),
|
||
Some(_) => return Err(Error::NotSupported(format!("orthotropic thermal expansion (material `{name}`) in explicit dynamics"))),
|
||
};
|
||
Ok(Material {
|
||
name: name.to_string(),
|
||
expansion,
|
||
elasticity: d,
|
||
density,
|
||
wave_speed: (stiffest / density).sqrt(),
|
||
modulus: modulus_bound(elastic, &d),
|
||
stiffest,
|
||
j2,
|
||
erosion: m.erosion.map(|e| e.plastic_strain),
|
||
})
|
||
}
|
||
|
||
/// The section's elements of a kind the explicit solver has, or why not.
|
||
fn check_formulation(kind: ElementKind, section_name: &str, f: &fieldcraft_model::Formulation) -> Result<()> {
|
||
let refuse = |why: &str| Err(Error::NotSupported(format!("section `{section_name}`: {why}")));
|
||
if f.hybrid {
|
||
return refuse("hybrid elements in explicit dynamics");
|
||
}
|
||
match kind {
|
||
ElementKind::Hex8 if f.incompatible_modes || f.integration == Integration::Full => refuse(
|
||
"fully integrated or incompatible-mode hex8 in explicit dynamics (explicit hex8 elements are one-point with hourglass control: set the section's integration to reduced)",
|
||
),
|
||
ElementKind::Hex8 => Ok(()),
|
||
ElementKind::Tet4 if f.integration == Integration::Reduced => refuse("reduced integration of tet4 elements"),
|
||
ElementKind::Tet4 => Ok(()),
|
||
ElementKind::Rod2 | ElementKind::Rod3 => Ok(()),
|
||
other => refuse(&format!("{} elements in explicit dynamics (it has one-point hex8 and nodal-pressure tet4 so far)", other.name())),
|
||
}
|
||
}
|
||
|
||
impl System {
|
||
pub fn new(model: &Model, resolved: &Resolved, step: &Step, settings: &ExplicitDynamic, max_dofs: usize) -> Result<System> {
|
||
if let Some(name) = model.couplings().keys().next() {
|
||
return Err(Error::NotSupported(format!("couplings (coupling `{name}`) in explicit dynamics")));
|
||
}
|
||
// Elements by label, with their section's material.
|
||
let mut materials: Vec<Material> = Vec::new();
|
||
let mut chosen: Vec<(u32, ElementKind, Vec<u32>, usize, Option<f64>)> = Vec::new();
|
||
for (name, section) in model.sections() {
|
||
let area = match section.kind {
|
||
SectionKind::Solid {} => None,
|
||
SectionKind::Rod { area } => Some(area),
|
||
_ => {
|
||
return Err(Error::NotSupported(format!(
|
||
"{} sections (section `{name}`); explicit dynamics has solid and truss elements so far",
|
||
section.kind.name()
|
||
)));
|
||
}
|
||
};
|
||
let mi = match materials.iter().position(|m| m.name == section.material) {
|
||
Some(i) => i,
|
||
None => {
|
||
materials.push(material(model, §ion.material)?);
|
||
materials.len() - 1
|
||
}
|
||
};
|
||
let set = resolved.sets.get(§ion.set).ok_or_else(|| invalid(format!("set `{}` is not resolved", section.set)))?;
|
||
for &label in &set.elements {
|
||
let el = model.element(label).ok_or_else(|| invalid(format!("element {label} is missing")))?;
|
||
if area.is_some() != matches!(el.kind, ElementKind::Rod2 | ElementKind::Rod3) {
|
||
return Err(invalid(format!(
|
||
"section `{name}`: a {} section cannot hold {} element {label}",
|
||
section.kind.name(),
|
||
el.kind.name()
|
||
)));
|
||
}
|
||
check_formulation(el.kind, name, §ion.formulation)?;
|
||
if el.kind == ElementKind::Tet4 && materials.get(mi).is_some_and(|m| m.j2.is_none()) {
|
||
return Err(Error::NotSupported(format!(
|
||
"nodal-pressure tet4 elements of the orthotropic material `{}` (section `{name}`)",
|
||
section.material
|
||
)));
|
||
}
|
||
chosen.push((label, el.kind, el.nodes.clone(), mi, area));
|
||
}
|
||
}
|
||
chosen.sort_by_key(|e| e.0);
|
||
let mut labels: Vec<u32> = chosen.iter().flat_map(|e| e.2.iter().copied()).collect();
|
||
// Mesh nodes only discrete elements use (reference points are refused by Discrete::new).
|
||
labels.extend(resolved.discrete_labels().filter(|l| model.node(*l).is_some()));
|
||
labels.sort_unstable();
|
||
labels.dedup();
|
||
let layout = DofLayout::uniform(labels.iter().copied(), DofSet::TRANSLATIONS);
|
||
if layout.dofs() > max_dofs {
|
||
return Err(Error::Limit(format!("the model has {} degrees of freedom; the limit is {max_dofs}", layout.dofs())));
|
||
}
|
||
if layout.stride() != Some(NODE_DOFS) {
|
||
return Err(Error::Internal("explicit dynamics needs three translations per node".into()));
|
||
}
|
||
if labels.is_empty() {
|
||
return Err(invalid("the model has no elements or point masses to analyse"));
|
||
}
|
||
let xyz: Vec<[f64; 3]> =
|
||
labels.iter().map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| invalid(format!("node {l} is missing")))).collect::<Result<_>>()?;
|
||
let (n, dofs) = (labels.len(), layout.dofs());
|
||
let mut sys = System {
|
||
layout,
|
||
elements: chosen.iter().map(|e| e.0).collect(),
|
||
materials,
|
||
hexes: Vec::new(),
|
||
tets: Vec::new(),
|
||
rods: Vec::new(),
|
||
mass: vec![0.0; n],
|
||
discrete: Discrete::default(),
|
||
discrete_increment: None,
|
||
prescribed: vec![None; dofs],
|
||
walls: Vec::new(),
|
||
initial_velocity: vec![0.0; dofs],
|
||
external: vec![0.0; dofs],
|
||
history_nodes: Vec::new(),
|
||
hourglass: settings.hourglass.kind,
|
||
hourglass_coefficient: settings.hourglass.coefficient,
|
||
nonlinear: settings.nonlinear_geometry,
|
||
stress_rate: settings.stress_rate,
|
||
warnings: resolved.warnings.iter().map(ToString::to_string).collect(),
|
||
labels,
|
||
xyz,
|
||
};
|
||
for (index, (label, kind, nodes, mi, area)) in chosen.into_iter().enumerate() {
|
||
let idx: Vec<usize> = nodes.iter().map(|l| sys.node(*l).ok_or_else(|| Error::Internal(format!("node {l}")))).collect::<Result<_>>()?;
|
||
let mat = sys.materials.get(mi).ok_or_else(|| Error::Internal("material index".into()))?;
|
||
let (rho, c) = (mat.density, mat.wave_speed);
|
||
let at = |e| Error::element(label, e);
|
||
match kind {
|
||
ElementKind::Hex8 => {
|
||
let nodes: [usize; 8] = idx.try_into().map_err(|_| invalid(format!("element {label} does not have 8 nodes")))?;
|
||
let x = nodes.map(|i| sys.xyz.get(i).copied().unwrap_or_default());
|
||
let k = Hex8::hex8(&x).map_err(at)?;
|
||
let gamma = k.hourglass_vectors(&x);
|
||
let mass = k.nodal_volume.map(|v| rho * v);
|
||
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 });
|
||
}
|
||
ElementKind::Tet4 => {
|
||
let nodes: [usize; 4] = idx.try_into().map_err(|_| invalid(format!("element {label} does not have 4 nodes")))?;
|
||
let x = nodes.map(|i| sys.xyz.get(i).copied().unwrap_or_default());
|
||
let k = Tet4::tet4(&x).map_err(at)?;
|
||
let mass = k.nodal_volume.map(|v| rho * v);
|
||
let omega2 = frequency_bound(&k, k.gradient_norm(), mat.modulus, &mass);
|
||
sys.lump(&nodes, &mass);
|
||
sys.tets.push(Tet {
|
||
label,
|
||
index,
|
||
nodes,
|
||
material: mi,
|
||
kinematics: k,
|
||
mass,
|
||
omega2,
|
||
length: 2.0 * c / omega2.sqrt(),
|
||
thermal_rate: 0.0,
|
||
});
|
||
}
|
||
ElementKind::Rod2 | ElementKind::Rod3 => {
|
||
let area = area.ok_or_else(|| Error::Internal("truss area".into()))?;
|
||
let source = model.material(&mat.name).ok_or_else(|| Error::Internal("truss material".into()))?;
|
||
if source.plastic.is_some() || source.damage.is_some() {
|
||
return Err(Error::NotSupported(format!(
|
||
"plastic or damaging materials on trusses (material `{}`, element {label}): explicit trusses are elastic so far",
|
||
mat.name
|
||
)));
|
||
}
|
||
if kind == ElementKind::Rod3 && settings.nonlinear_geometry {
|
||
return Err(Error::NotSupported(format!(
|
||
"three-node trusses with nonlinear geometry (element {label}); set the step's nonlinear_geometry to false"
|
||
)));
|
||
}
|
||
let young =
|
||
fieldcraft_fem::truss::young_modulus(source.elastic.as_ref().ok_or_else(|| Error::Internal("truss material".into()))?)
|
||
.map_err(|e| invalid(format!("material `{}` of element {label}: {e}", mat.name)))?;
|
||
let x: Vec<[f64; 3]> = idx.iter().map(|i| sys.xyz.get(*i).copied().unwrap_or_default()).collect();
|
||
let rod = Rod::new(label, index, kind, idx.clone(), &x, mi, crate::discrete::RodSection { young, density: rho, area })?;
|
||
let mass: Vec<f64> = rod.nodal_volume.iter().map(|v| rho * v).collect();
|
||
sys.lump(&idx, &mass);
|
||
sys.rods.push(rod);
|
||
}
|
||
other => return Err(Error::Internal(format!("{} element kind", other.name()))),
|
||
}
|
||
}
|
||
sys.discrete = Discrete::new(resolved, sys.dofs(), &sys.xyz, |l| sys.node(l))?;
|
||
sys.discrete_increment = sys.discrete.critical_increment(&sys.dof_mass(), sys.nonlinear);
|
||
sys.prescribe(model, resolved, step, settings.end_time)?;
|
||
sys.load(model, resolved, step)?;
|
||
sys.temperatures(model, resolved, step, settings.end_time)?;
|
||
sys.conditions(resolved, settings)?;
|
||
Ok(sys)
|
||
}
|
||
|
||
/// Add an element's row-sum lumped masses.
|
||
fn lump(&mut self, nodes: &[usize], mass: &[f64]) {
|
||
for (node, v) in nodes.iter().zip(mass) {
|
||
if let Some(m) = self.mass.get_mut(*node) {
|
||
*m += v;
|
||
}
|
||
}
|
||
}
|
||
|
||
pub fn node(&self, label: u32) -> Option<usize> {
|
||
self.labels.binary_search(&label).ok()
|
||
}
|
||
|
||
pub fn dofs(&self) -> usize {
|
||
self.layout.dofs()
|
||
}
|
||
|
||
/// The lumped mass per degree of freedom: the elements' and the point masses'.
|
||
pub fn dof_mass(&self) -> Vec<f64> {
|
||
(0..self.dofs())
|
||
.map(|d| self.mass.get(crate::discrete::node_direction(d).0).copied().unwrap_or(0.0) + self.discrete.mass.get(d).copied().unwrap_or(0.0))
|
||
.collect()
|
||
}
|
||
|
||
/// Prescribed velocities from the step's boundary conditions: zero for held degrees of
|
||
/// freedom, value / end time for a displacement reached over the step.
|
||
fn prescribe(&mut self, model: &Model, resolved: &Resolved, step: &Step, end_time: f64) -> Result<()> {
|
||
for name in &step.boundary_conditions {
|
||
let bc = model.boundary_condition(name).ok_or_else(|| invalid(format!("boundary condition `{name}` is missing")))?;
|
||
if let Target::ReferencePoint(rp) = bc.target() {
|
||
return Err(Error::NotSupported(format!("boundary conditions on reference points (`{name}` on `{rp}`)")));
|
||
}
|
||
if bc.is_thermal() {
|
||
return Err(Error::NotSupported(format!("temperature boundary conditions (`{name}`) in explicit dynamics")));
|
||
}
|
||
for &label in resolved.target_nodes(bc.target()) {
|
||
let Some(node) = self.node(label) else { continue };
|
||
for kind in self.layout.kinds(node).iter() {
|
||
let dof = kind.number();
|
||
let Some(value) = bc.value(dof) else { continue };
|
||
let velocity = value / end_time;
|
||
if !velocity.is_finite() {
|
||
return Err(invalid(format!("boundary condition `{name}`: the velocity {value} / {end_time} is not finite")));
|
||
}
|
||
let index = self.layout.dof(node, kind).ok_or_else(|| Error::Internal("dof index".into()))?;
|
||
let slot = self.prescribed.get_mut(index).ok_or_else(|| Error::Internal("dof index".into()))?;
|
||
match *slot {
|
||
Some(old) if old != velocity => {
|
||
return Err(invalid(format!(
|
||
"node {label} direction {dof} moves at both {old} and {velocity} (boundary condition `{name}` and another)"
|
||
)));
|
||
}
|
||
_ => *slot = Some(velocity),
|
||
}
|
||
}
|
||
}
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
/// The thermal strain rates of the step's temperature loads: each element's mean nodal
|
||
/// temperature change from its material's reference (nodes no load reaches stay at it) is
|
||
/// reached linearly over the step, like a prescribed displacement, so the thermal strain
|
||
/// α ΔT grows at the rate α ΔT / end time (no thermal shock at the start).
|
||
fn temperatures(&mut self, model: &Model, resolved: &Resolved, step: &Step, end_time: f64) -> Result<()> {
|
||
let mut field: Vec<Option<f64>> = vec![None; self.labels.len()];
|
||
let mut any = false;
|
||
for name in &step.loads {
|
||
let Some(Load::Temperature { target, value }) = model.load(name) else { continue };
|
||
any = true;
|
||
for &label in resolved.target_nodes(target) {
|
||
let Some(slot) = self.node(label).and_then(|n| field.get_mut(n)) else { continue };
|
||
match *slot {
|
||
Some(old) if old != *value => {
|
||
return Err(invalid(format!("node {label} is given both {old} and {value} (temperature load `{name}` and another)")));
|
||
}
|
||
_ => *slot = Some(*value),
|
||
}
|
||
}
|
||
}
|
||
if !any {
|
||
return Ok(());
|
||
}
|
||
if !self.rods.is_empty() {
|
||
self.warnings
|
||
.push("trusses take no thermal strain in explicit dynamics yet; the temperature loads strain the solid elements only".into());
|
||
}
|
||
let materials = &self.materials;
|
||
let rate = |nodes: &[usize], material: usize| -> f64 {
|
||
let Some((alpha, reference)) = materials.get(material).and_then(|m| m.expansion) else { return 0.0 };
|
||
let dt: f64 = nodes.iter().map(|n| field.get(*n).copied().flatten().map_or(0.0, |t| t - reference)).sum::<f64>() / nodes.len() as f64;
|
||
alpha * dt / end_time
|
||
};
|
||
for h in &mut self.hexes {
|
||
h.thermal_rate = rate(&h.nodes, h.material);
|
||
}
|
||
for t in &mut self.tets {
|
||
t.thermal_rate = rate(&t.nodes, t.material);
|
||
}
|
||
if self.hexes.iter().map(|h| h.thermal_rate).chain(self.tets.iter().map(|t| t.thermal_rate)).any(|r| !r.is_finite()) {
|
||
return Err(invalid("the thermal strain rate is not finite"));
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
/// The constant load vector: forces and tractions as the model distributes them, consistent
|
||
/// pressure loads, gravity and body forces with the lumped (row-sum) weights.
|
||
fn load(&mut self, model: &Model, resolved: &Resolved, step: &Step) -> Result<()> {
|
||
let mut f = std::mem::take(&mut self.external);
|
||
for name in &step.loads {
|
||
let load = model.load(name).ok_or_else(|| invalid(format!("load `{name}` is missing")))?;
|
||
let off_mesh = |l: u32| invalid(format!("load `{name}`: node {l} belongs to no element"));
|
||
match load {
|
||
Load::Force { target: Target::ReferencePoint(rp), .. } => {
|
||
return Err(Error::NotSupported(format!("loads on reference points (`{name}` on `{rp}`)")));
|
||
}
|
||
Load::Moment { .. } => {
|
||
return Err(Error::NotSupported(format!("moments (`{name}`): solid nodes have no rotations")));
|
||
}
|
||
Load::SurfaceFlux { .. } | Load::Convection { .. } | Load::HeatSource { .. } => {
|
||
return Err(Error::NotSupported(format!("{} loads (`{name}`) in explicit dynamics", load.kind())));
|
||
}
|
||
// Thermal strains, set by `temperatures`.
|
||
Load::Temperature { .. } => {}
|
||
Load::Force { .. } | Load::Traction { .. } => {
|
||
for (label, v) in resolved.nodal_loads.get(name).map(Vec::as_slice).unwrap_or_default() {
|
||
let node = self.node(*label).ok_or_else(|| off_mesh(*label))?;
|
||
add_at(&mut f, node, v);
|
||
}
|
||
}
|
||
Load::Pressure { set, magnitude } => {
|
||
for fc in resolved.sets.get(set).map(|s| s.faces.as_slice()).unwrap_or_default() {
|
||
let labels = model
|
||
.face_nodes(*fc)
|
||
.ok_or_else(|| invalid(format!("load `{name}`: element {} face {} is missing", fc.element, fc.face)))?;
|
||
let idx: Vec<usize> = labels.iter().map(|l| self.node(*l).ok_or_else(|| off_mesh(*l))).collect::<Result<_>>()?;
|
||
let xyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| self.xyz.get(*i).copied()).collect();
|
||
let forces = face::pressure(&xyz, *magnitude).map_err(|e| Error::element(fc.element, e))?;
|
||
for (i, v) in idx.iter().zip(forces) {
|
||
add_at(&mut f, *i, &v);
|
||
}
|
||
}
|
||
}
|
||
Load::Gravity { acceleration, set } => {
|
||
let chosen: Option<&[u32]> = set.as_ref().map(|s| resolved.sets.get(s).map(|r| r.elements.as_slice()).unwrap_or_default());
|
||
self.each_element(chosen, |nodes, weights, rho| {
|
||
for (node, w) in nodes.iter().zip(weights) {
|
||
add_at(&mut f, *node, &acceleration.map(|a| rho * w * a));
|
||
}
|
||
});
|
||
// Point masses take gravity on the whole model.
|
||
if set.is_none() {
|
||
for (d, (x, m)) in f.iter_mut().zip(&self.discrete.mass).enumerate() {
|
||
*x += m * acceleration.get(crate::discrete::node_direction(d).1).copied().unwrap_or(0.0);
|
||
}
|
||
}
|
||
}
|
||
Load::BodyForce { set, vector } => {
|
||
let chosen = resolved.sets.get(set).map(|r| r.elements.as_slice()).unwrap_or_default();
|
||
self.each_element(Some(chosen), |nodes, weights, _| {
|
||
for (node, w) in nodes.iter().zip(weights) {
|
||
add_at(&mut f, *node, &vector.map(|b| w * b));
|
||
}
|
||
});
|
||
}
|
||
Load::Centrifugal { angular_velocity, point, axis, set } => {
|
||
let chosen: Option<&[u32]> = set.as_ref().map(|s| resolved.sets.get(s).map(|r| r.elements.as_slice()).unwrap_or_default());
|
||
let xyz = &self.xyz;
|
||
self.each_element(chosen, |nodes, weights, rho| {
|
||
for (node, w) in nodes.iter().zip(weights) {
|
||
let x = xyz.get(*node).copied().unwrap_or_default();
|
||
let a = fieldcraft_model::centrifugal_acceleration(*angular_velocity, *point, *axis, x);
|
||
add_at(&mut f, *node, &a.map(|c| rho * w * c));
|
||
}
|
||
});
|
||
}
|
||
}
|
||
}
|
||
self.external = f;
|
||
Ok(())
|
||
}
|
||
|
||
/// Visit the elements (all, or those listed): their nodes, nodal volumes and density.
|
||
fn each_element(&self, only: Option<&[u32]>, mut visit: impl FnMut(&[usize], &[f64], f64)) {
|
||
let wanted = |label: u32| only.is_none_or(|l| l.binary_search(&label).is_ok());
|
||
let density = |m: usize| self.materials.get(m).map_or(0.0, |m| m.density);
|
||
for h in self.hexes.iter().filter(|h| wanted(h.label)) {
|
||
visit(&h.nodes, &h.kinematics.nodal_volume, density(h.material));
|
||
}
|
||
for t in self.tets.iter().filter(|t| wanted(t.label)) {
|
||
visit(&t.nodes, &t.kinematics.nodal_volume, density(t.material));
|
||
}
|
||
for r in self.rods.iter().filter(|r| wanted(r.label)) {
|
||
visit(&r.nodes, &r.nodal_volume, density(r.material));
|
||
}
|
||
}
|
||
|
||
/// The analysed nodes of a resolved set.
|
||
fn nodes_of(&self, resolved: &Resolved, set: &str) -> Vec<usize> {
|
||
resolved.sets.get(set).map(|s| s.nodes.as_slice()).unwrap_or_default().iter().filter_map(|l| self.node(*l)).collect()
|
||
}
|
||
|
||
/// Initial velocities, rigid walls and history nodes of the step.
|
||
fn conditions(&mut self, resolved: &Resolved, settings: &ExplicitDynamic) -> Result<()> {
|
||
for iv in &settings.initial_velocities {
|
||
for node in self.nodes_of(resolved, &iv.set) {
|
||
for (k, c) in iv.velocity.iter().enumerate() {
|
||
if let Some(v) =
|
||
DofKind::translation(k).and_then(|kind| self.layout.dof(node, kind)).and_then(|d| self.initial_velocity.get_mut(d))
|
||
{
|
||
*v = *c;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
// A prescribed velocity holds from the start.
|
||
for (v, p) in self.initial_velocity.iter_mut().zip(&self.prescribed) {
|
||
if let Some(p) = p {
|
||
*v = *p;
|
||
}
|
||
}
|
||
let size = self.xyz.iter().flatten().fold(0.0f64, |m, c| m.max(c.abs())).max(1.0);
|
||
for w in &settings.rigid_walls {
|
||
let len = w.normal.iter().map(|c| c * c).sum::<f64>().sqrt();
|
||
let normal = w.normal.map(|c| c / len);
|
||
let nodes = self.nodes_of(resolved, &w.set);
|
||
for &i in &nodes {
|
||
let x = self.xyz.get(i).copied().unwrap_or_default();
|
||
let gap: f64 = (0..3).map(|k| (x[k] - w.point[k]) * normal[k]).sum();
|
||
if gap < -1e-9 * size {
|
||
let label = self.labels.get(i).copied().unwrap_or(0);
|
||
return Err(invalid(format!("node {label} of set `{}` starts behind its rigid wall (gap {gap:e})", w.set)));
|
||
}
|
||
}
|
||
self.walls.push(Wall { point: w.point, normal, nodes });
|
||
}
|
||
if let Some(set) = &settings.history_nodes {
|
||
let nodes = self.nodes_of(resolved, set);
|
||
if nodes.len() > MAX_HISTORY_NODES {
|
||
return Err(Error::Limit(format!("set `{set}` has {} nodes; histories are recorded for at most {MAX_HISTORY_NODES}", nodes.len())));
|
||
}
|
||
self.history_nodes = nodes;
|
||
}
|
||
Ok(())
|
||
}
|
||
}
|