Compare commits
2
Commits
cca29aac8f
...
e30cfe4ce9
| Author | SHA1 | Date | |
|---|---|---|---|
 Omar SobhandClaude Opus 5
|
e30cfe4ce9 | ||
|
Omar SobhandClaude Opus 5
|
4c2cea36aa |
@@ -1,13 +1,25 @@
|
||||
// Copyright (c) 2024 RustyTorch++ Team
|
||||
// Licensed under the Apache License, Version 2.0
|
||||
|
||||
//! Modal analysis for eigenvalue problems.
|
||||
//! Modal analysis: natural frequencies and mode shapes.
|
||||
//!
|
||||
//! Solves the generalized eigenvalue problem `K φ = λ M φ` on the free
|
||||
//! degrees of freedom, where `λ = ω²`.
|
||||
//!
|
||||
//! This previously returned `DVector::zeros(num_modes)` without assembling
|
||||
//! anything. Both halves it needs already existed —
|
||||
//! [`crate::assembly::GlobalAssembler`] for `K` and `M`, and
|
||||
//! [`crate::solvers::eigenvalue::EigenvalueSolver`] for the eigenproblem —
|
||||
//! they were simply never connected.
|
||||
|
||||
use super::{Analysis, AnalysisConfig, AnalysisData, AnalysisResults};
|
||||
use crate::error::FeaResult;
|
||||
use crate::assembly::{AdvancedDofNumbering, DofMappingStrategy, GlobalAssembler};
|
||||
use crate::boundary::{BoundaryCondition, BoundaryConditionSet};
|
||||
use crate::error::{AnalysisError, FeaResult};
|
||||
use crate::materials::MaterialDatabase;
|
||||
use crate::mesh::Mesh;
|
||||
use nalgebra::DVector;
|
||||
use crate::solvers::eigenvalue::{EigenvalueSolver, ModalResults};
|
||||
use nalgebra::{DMatrix, DVector};
|
||||
|
||||
/// Modal analysis for computing natural frequencies and mode shapes.
|
||||
#[derive(Debug)]
|
||||
@@ -16,6 +28,8 @@ pub struct ModalAnalysis {
|
||||
materials: MaterialDatabase,
|
||||
num_modes: usize,
|
||||
config: AnalysisConfig,
|
||||
boundary_conditions: BoundaryConditionSet,
|
||||
shift: Option<f64>,
|
||||
progress: f64,
|
||||
complete: bool,
|
||||
}
|
||||
@@ -32,27 +46,153 @@ impl ModalAnalysis {
|
||||
materials,
|
||||
num_modes,
|
||||
config,
|
||||
boundary_conditions: BoundaryConditionSet::new(),
|
||||
shift: None,
|
||||
progress: 0.0,
|
||||
complete: false,
|
||||
}
|
||||
}
|
||||
|
||||
/// Constrain the structure.
|
||||
///
|
||||
/// Without this the stiffness matrix is singular: an unconstrained body
|
||||
/// has rigid-body modes at zero frequency. That is a physical fact rather
|
||||
/// than a numerical accident, so [`Self::run`] reports it as an error
|
||||
/// naming the remedy instead of returning six meaningless near-zero
|
||||
/// modes.
|
||||
pub fn with_boundary_conditions(mut self, boundary_conditions: BoundaryConditionSet) -> Self {
|
||||
self.boundary_conditions = boundary_conditions;
|
||||
self
|
||||
}
|
||||
|
||||
/// Solve for the modes nearest this frequency-squared value instead of
|
||||
/// the lowest ones.
|
||||
///
|
||||
/// The one case that genuinely needs it is a free-free structure, where a
|
||||
/// small negative shift moves the rigid-body modes off the singularity.
|
||||
pub fn with_shift(mut self, shift: f64) -> Self {
|
||||
self.shift = Some(shift);
|
||||
self
|
||||
}
|
||||
|
||||
/// Number the degrees of freedom and apply displacement constraints.
|
||||
fn setup_dof_numbering(&self) -> FeaResult<AdvancedDofNumbering> {
|
||||
let mut dof_numbering = AdvancedDofNumbering::displacement_only(
|
||||
&self.mesh,
|
||||
DofMappingStrategy::BandwidthOptimized,
|
||||
)?;
|
||||
|
||||
for bc in self.boundary_conditions.conditions() {
|
||||
// Only Dirichlet conditions remove degrees of freedom. A modal
|
||||
// analysis has no applied load, so Neumann conditions have
|
||||
// nothing to contribute and are ignored rather than silently
|
||||
// mis-applied.
|
||||
if let BoundaryCondition::Dirichlet(dirichlet) = bc {
|
||||
for &node_id in &dirichlet.nodes {
|
||||
for &component in &dirichlet.components {
|
||||
if let Some(dof) = dof_numbering.get_dof(node_id, component) {
|
||||
dof_numbering.constrain_dof(dof)?;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Ok(dof_numbering)
|
||||
}
|
||||
}
|
||||
|
||||
impl Analysis for ModalAnalysis {
|
||||
fn run(&mut self) -> FeaResult<AnalysisResults> {
|
||||
// Simplified modal analysis
|
||||
let num_dofs = self.mesh.num_nodes() * 3;
|
||||
let solution = DVector::zeros(num_dofs);
|
||||
if self.num_modes == 0 {
|
||||
return Err(AnalysisError::InvalidConfiguration(
|
||||
"modal analysis requires at least one mode".to_string(),
|
||||
)
|
||||
.into());
|
||||
}
|
||||
if self.mesh.num_elements() == 0 {
|
||||
return Err(AnalysisError::InvalidConfiguration(
|
||||
"mesh contains no elements".to_string(),
|
||||
)
|
||||
.into());
|
||||
}
|
||||
|
||||
let dof_numbering = self.setup_dof_numbering()?;
|
||||
self.progress = 0.2;
|
||||
|
||||
let assembler =
|
||||
GlobalAssembler::new(self.config.assembly_options.clone(), self.materials.clone());
|
||||
let system =
|
||||
assembler.assemble_system(&self.mesh, &dof_numbering.to_dof_numbering(), 0.0)?;
|
||||
self.progress = 0.5;
|
||||
|
||||
let mass = system.mass_matrix.as_ref().ok_or_else(|| {
|
||||
AnalysisError::InvalidConfiguration(
|
||||
"assembly produced no mass matrix; modal analysis needs one".to_string(),
|
||||
)
|
||||
})?;
|
||||
|
||||
let free_dofs = &system.dof_numbering.free_dofs;
|
||||
if free_dofs.is_empty() {
|
||||
return Err(AnalysisError::InvalidConfiguration(
|
||||
"every degree of freedom is constrained; there is nothing to vibrate".to_string(),
|
||||
)
|
||||
.into());
|
||||
}
|
||||
if self.num_modes > free_dofs.len() {
|
||||
return Err(AnalysisError::InvalidConfiguration(format!(
|
||||
"requested {} modes from a system with {} free degrees of freedom",
|
||||
self.num_modes,
|
||||
free_dofs.len()
|
||||
))
|
||||
.into());
|
||||
}
|
||||
|
||||
let free_stiffness = system
|
||||
.stiffness_matrix
|
||||
.extract_submatrix(free_dofs, free_dofs)?;
|
||||
let free_mass = mass.extract_submatrix(free_dofs, free_dofs)?;
|
||||
self.progress = 0.7;
|
||||
|
||||
let mut solver = EigenvalueSolver::new(self.num_modes);
|
||||
if let Some(shift) = self.shift {
|
||||
solver = solver.with_shift(shift);
|
||||
}
|
||||
let (eigenvalues, free_shapes) = solver.solve(&free_stiffness, &free_mass)?;
|
||||
self.progress = 0.9;
|
||||
|
||||
// Expand each mode shape from the free-DOF space back to the full
|
||||
// DOF vector, leaving constrained entries at zero.
|
||||
let total_dofs = system.dof_numbering.total_dofs;
|
||||
let mut mode_shapes = DMatrix::zeros(total_dofs, self.num_modes);
|
||||
for mode in 0..self.num_modes {
|
||||
for (i, &global_dof) in free_dofs.iter().enumerate() {
|
||||
mode_shapes[(global_dof, mode)] = free_shapes[(i, mode)];
|
||||
}
|
||||
}
|
||||
|
||||
let modal = ModalResults::from_eigenvalues(eigenvalues, free_shapes);
|
||||
|
||||
// The reported solution is the fundamental mode shape, so that a
|
||||
// caller treating this like any other analysis gets something
|
||||
// physical rather than a zero vector.
|
||||
let solution =
|
||||
DVector::from_iterator(total_dofs, (0..total_dofs).map(|d| mode_shapes[(d, 0)]));
|
||||
|
||||
let mut results = AnalysisResults::new("Modal".to_string(), solution);
|
||||
|
||||
// Add modal data
|
||||
let eigenvalues = DVector::zeros(self.num_modes);
|
||||
let frequencies: DVector<f64> =
|
||||
eigenvalues.map(|lambda: f64| (lambda / (2.0 * std::f64::consts::PI)).sqrt());
|
||||
|
||||
results.add_data("eigenvalues".to_string(), AnalysisData::Vector(eigenvalues));
|
||||
results.add_data("frequencies".to_string(), AnalysisData::Vector(frequencies));
|
||||
results.add_data(
|
||||
"eigenvalues".to_string(),
|
||||
AnalysisData::Vector(modal.eigenvalues.clone()),
|
||||
);
|
||||
results.add_data(
|
||||
"frequencies".to_string(),
|
||||
AnalysisData::Vector(modal.frequencies.clone()),
|
||||
);
|
||||
results.add_data(
|
||||
"periods".to_string(),
|
||||
AnalysisData::Vector(modal.periods.clone()),
|
||||
);
|
||||
results.add_data("mode_shapes".to_string(), AnalysisData::Matrix(mode_shapes));
|
||||
|
||||
self.progress = 1.0;
|
||||
self.complete = true;
|
||||
|
||||
@@ -21,7 +21,12 @@ pub enum DofMappingStrategy {
|
||||
}
|
||||
|
||||
/// DOF component types for different physics.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
|
||||
///
|
||||
/// The declaration order is the **canonical DOF order within a node** — the
|
||||
/// order element matrices are built in, so `[u, v, w]` before rotations
|
||||
/// before scalar fields. [`DofComponent::canonical_index`] exposes it, and
|
||||
/// anything that collects a node's degrees of freedom must sort by it.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
|
||||
pub enum DofComponent {
|
||||
/// Displacement in X direction
|
||||
DisplacementX,
|
||||
@@ -46,6 +51,26 @@ pub enum DofComponent {
|
||||
}
|
||||
|
||||
impl DofComponent {
|
||||
/// Position of this component in the canonical within-node DOF order.
|
||||
///
|
||||
/// Element matrices interleave displacement components per node —
|
||||
/// `[u₀, v₀, u₁, v₁, …]` — so any mapping from a node to its global DOFs
|
||||
/// must present them in this order or the assembly is transposed.
|
||||
pub fn canonical_index(self) -> usize {
|
||||
match self {
|
||||
Self::DisplacementX => 0,
|
||||
Self::DisplacementY => 1,
|
||||
Self::DisplacementZ => 2,
|
||||
Self::RotationX => 3,
|
||||
Self::RotationY => 4,
|
||||
Self::RotationZ => 5,
|
||||
Self::Temperature => 6,
|
||||
Self::Pressure => 7,
|
||||
Self::ElectricPotential => 8,
|
||||
Self::MagneticPotential => 9,
|
||||
}
|
||||
}
|
||||
|
||||
/// Get all displacement components.
|
||||
pub fn displacement_components() -> Vec<Self> {
|
||||
vec![
|
||||
@@ -115,10 +140,22 @@ impl AdvancedDofNumbering {
|
||||
}
|
||||
|
||||
/// Create simple displacement-only DOF numbering.
|
||||
///
|
||||
/// The number of displacement components per node follows the **mesh's
|
||||
/// spatial dimension**: two for a planar mesh, three for a solid one.
|
||||
///
|
||||
/// This previously always allocated all three components. On a 2-D mesh
|
||||
/// that gave a node three degrees of freedom while the element matrices
|
||||
/// — sized from `ElementType::spatial_dimension()` — supplied two, so
|
||||
/// assembly rejected every contribution with a dimension mismatch. A
|
||||
/// planar mesh has no out-of-plane displacement to number, and numbering
|
||||
/// one produces a global stiffness matrix with an entirely zero row and
|
||||
/// column, which is singular by construction.
|
||||
pub fn displacement_only(mesh: &Mesh, strategy: DofMappingStrategy) -> FeaResult<Self> {
|
||||
let components = DofComponent::displacement_components();
|
||||
let mut node_components = HashMap::new();
|
||||
let components = components[..mesh.spatial_dimension.min(components.len())].to_vec();
|
||||
|
||||
let mut node_components = HashMap::new();
|
||||
for &node_id in mesh.nodes.keys() {
|
||||
node_components.insert(node_id, components.clone());
|
||||
}
|
||||
@@ -473,16 +510,46 @@ impl AdvancedDofNumbering {
|
||||
}
|
||||
|
||||
/// Convert to simplified DofNumbering for compatibility with legacy code.
|
||||
/// Flatten into the simple [`super::DofNumbering`] the assembler consumes.
|
||||
///
|
||||
/// Each node's degrees of freedom are ordered by
|
||||
/// [`DofComponent`]'s canonical order — `[u, v, w]`, then rotations —
|
||||
/// because `GlobalAssembler` maps element local DOF `k` onto
|
||||
/// `get_node_dofs(node)[k]`, and element matrices are built in exactly
|
||||
/// that order.
|
||||
///
|
||||
/// This previously pushed the DOFs in `HashMap` iteration order, which is
|
||||
/// unspecified. When it came out as `[v, u]` the assembler wrote the
|
||||
/// element's `u` row into the global `v` row: the resulting matrix was
|
||||
/// still symmetric, still had the right rigid-body null space and still
|
||||
/// summed to the right total mass, but described a structure whose axes
|
||||
/// were transposed per node — and `get_dof(node, DisplacementX)` then
|
||||
/// pointed at the wrong row, so constraints were applied to the wrong
|
||||
/// direction too.
|
||||
pub fn to_dof_numbering(&self) -> super::DofNumbering {
|
||||
use super::DofNumbering;
|
||||
let mut node_to_dofs: HashMap<NodeId, Vec<usize>> = HashMap::new();
|
||||
let mut node_components: HashMap<NodeId, Vec<(DofComponent, usize)>> = HashMap::new();
|
||||
let mut dof_to_node: HashMap<usize, NodeId> = HashMap::new();
|
||||
|
||||
for ((node_id, _component), &dof) in &self.node_component_to_dof {
|
||||
node_to_dofs.entry(*node_id).or_default().push(dof);
|
||||
for ((node_id, component), &dof) in &self.node_component_to_dof {
|
||||
node_components
|
||||
.entry(*node_id)
|
||||
.or_default()
|
||||
.push((*component, dof));
|
||||
dof_to_node.insert(dof, *node_id);
|
||||
}
|
||||
|
||||
let node_to_dofs: HashMap<NodeId, Vec<usize>> = node_components
|
||||
.into_iter()
|
||||
.map(|(node_id, mut components)| {
|
||||
components.sort_by_key(|(component, _)| component.canonical_index());
|
||||
(
|
||||
node_id,
|
||||
components.into_iter().map(|(_, dof)| dof).collect(),
|
||||
)
|
||||
})
|
||||
.collect();
|
||||
|
||||
DofNumbering {
|
||||
total_dofs: self.total_dofs,
|
||||
node_to_dofs,
|
||||
|
||||
@@ -219,13 +219,40 @@ impl JacobianEval {
|
||||
}
|
||||
|
||||
/// Transform derivatives from natural to physical coordinates.
|
||||
/// Map shape function derivatives from natural to physical coordinates.
|
||||
///
|
||||
/// `natural_derivatives` is `(num_nodes × parametric_dim)` — the layout
|
||||
/// [`crate::elements::ShapeFunctionEval`] produces — and the result is
|
||||
/// `(num_nodes × spatial_dim)`, which is what
|
||||
/// [`ElementMatrixComputer::strain_displacement_matrix`] indexes as
|
||||
/// `derivatives[(node, direction)]`.
|
||||
///
|
||||
/// By the chain rule `∂Nᵢ/∂x_k = Σⱼ (∂Nᵢ/∂ξⱼ)(∂ξⱼ/∂x_k)`, and
|
||||
/// `J⁻¹[(j, k)] = ∂ξⱼ/∂x_k`, so the product is
|
||||
/// `natural_derivatives · J⁻¹` — not `J⁻ᵀ · natural_derivatives`, which
|
||||
/// this previously computed. The two agree only when both matrices are
|
||||
/// square and symmetric; for any element with more nodes than parametric
|
||||
/// directions — that is, every element — the old form was a dimension
|
||||
/// mismatch and panicked inside the BLAS call.
|
||||
pub fn transform_derivatives(
|
||||
&self,
|
||||
natural_derivatives: &DMatrix<f64>,
|
||||
) -> FeaResult<DMatrix<f64>> {
|
||||
// Physical derivatives = inverse_jacobian^T * natural_derivatives
|
||||
let physical_derivatives = self.inverse_jacobian.transpose() * natural_derivatives;
|
||||
Ok(physical_derivatives)
|
||||
if natural_derivatives.ncols() != self.inverse_jacobian.nrows() {
|
||||
return Err(ElementError::InvalidGeometry {
|
||||
message: format!(
|
||||
"shape derivative layout is ({}x{}), expected {} columns to \
|
||||
match the {}-dimensional parametric space",
|
||||
natural_derivatives.nrows(),
|
||||
natural_derivatives.ncols(),
|
||||
self.inverse_jacobian.nrows(),
|
||||
self.inverse_jacobian.nrows(),
|
||||
),
|
||||
}
|
||||
.into());
|
||||
}
|
||||
|
||||
Ok(natural_derivatives * &self.inverse_jacobian)
|
||||
}
|
||||
|
||||
/// Get the Jacobian determinant (used for integration).
|
||||
|
||||
@@ -321,25 +321,76 @@ impl StandardFiniteElement {
|
||||
}
|
||||
}
|
||||
|
||||
/// Compute element matrices (stiffness, mass, force)
|
||||
/// Compute the element stiffness, mass and force matrices.
|
||||
///
|
||||
/// This is the function [`crate::assembly::GlobalAssembler`] calls for
|
||||
/// every element in the mesh, so every global matrix in the crate is
|
||||
/// downstream of it. It integrates `K = ∫ Bᵀ D B dV` and
|
||||
/// `M = ∫ ρ N ᵀN dV` by the element's own quadrature rule, via
|
||||
/// [`ElementMatrixComputer`].
|
||||
///
|
||||
/// # Degree-of-freedom layout
|
||||
///
|
||||
/// Displacement components are interleaved per node — `[u₀, v₀, u₁, v₁,
|
||||
/// …]` in 2-D — matching the strain-displacement matrix assembled by
|
||||
/// [`ElementMatrixComputer::compute_stiffness_matrix`] and the node DOF
|
||||
/// blocks handed out by `DofNumbering`.
|
||||
///
|
||||
/// `compute_mass_matrix` integrates the *scalar* field mass,
|
||||
/// `num_nodes × num_nodes`, because `∫ ρ Nᵢ Nⱼ dV` does not depend on the
|
||||
/// component. The vector-valued mass matrix is its Kronecker product with
|
||||
/// the spatial identity: each displacement component carries the same
|
||||
/// mass and the components do not couple. Skipping that expansion is how
|
||||
/// a mass matrix ends up dimensionally inconsistent with its stiffness.
|
||||
///
|
||||
/// # Force vector
|
||||
///
|
||||
/// Zero. Body forces and surface tractions enter through the boundary
|
||||
/// condition set rather than here, so an element with no body load
|
||||
/// contributes no force — which is different from the previous behaviour
|
||||
/// of returning zero because nothing was computed at all.
|
||||
pub fn compute_element_matrices(
|
||||
&self,
|
||||
_material: &dyn crate::materials::Material,
|
||||
material: &dyn crate::materials::Material,
|
||||
_time: f64,
|
||||
) -> FeaResult<ElementMatrix> {
|
||||
let properties = material.properties();
|
||||
let num_nodes = self.num_nodes();
|
||||
let dofs_per_node = self.spatial_dimension();
|
||||
let total_dofs = num_nodes * dofs_per_node;
|
||||
let dim = self.spatial_dimension();
|
||||
let total_dofs = num_nodes * dim;
|
||||
|
||||
// Create placeholder matrices
|
||||
let stiffness_matrix = nalgebra::DMatrix::zeros(total_dofs, total_dofs);
|
||||
let force_vector = nalgebra::DVector::zeros(total_dofs);
|
||||
let mass_matrix = Some(nalgebra::DMatrix::zeros(total_dofs, total_dofs));
|
||||
let stiffness_matrix = ElementMatrixComputer::compute_stiffness_matrix(
|
||||
self,
|
||||
&self.node_coordinates,
|
||||
properties.elastic_modulus,
|
||||
properties.poisson_ratio,
|
||||
None,
|
||||
)?
|
||||
.matrix;
|
||||
|
||||
let scalar_mass = ElementMatrixComputer::compute_consistent_mass_matrix(
|
||||
self,
|
||||
&self.node_coordinates,
|
||||
properties.density,
|
||||
None,
|
||||
)?
|
||||
.matrix;
|
||||
|
||||
// M_vector = M_scalar ⊗ I_dim, interleaved to match the DOF layout.
|
||||
let mut mass_matrix = nalgebra::DMatrix::zeros(total_dofs, total_dofs);
|
||||
for i in 0..num_nodes {
|
||||
for j in 0..num_nodes {
|
||||
let m_ij = scalar_mass[(i, j)];
|
||||
for d in 0..dim {
|
||||
mass_matrix[(i * dim + d, j * dim + d)] = m_ij;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Ok(ElementMatrix {
|
||||
stiffness_matrix,
|
||||
force_vector,
|
||||
mass_matrix,
|
||||
force_vector: nalgebra::DVector::zeros(total_dofs),
|
||||
mass_matrix: Some(mass_matrix),
|
||||
})
|
||||
}
|
||||
|
||||
@@ -411,14 +462,72 @@ impl FiniteElement for StandardFiniteElement {
|
||||
jacobian::JacobianEval::compute(self, coords, node_coords)
|
||||
}
|
||||
|
||||
/// Select the quadrature rule for this element's reference domain.
|
||||
///
|
||||
/// This previously returned `QuadratureRule::new(vec![], ..)` — a rule
|
||||
/// with **no points**. Every integration loop in the crate iterates over
|
||||
/// `rule.points`, so an empty rule does not fail: it silently skips the
|
||||
/// loop body and yields a zero matrix. That is why the element stiffness
|
||||
/// and mass matrices were zero, and it is the reason an unsupported
|
||||
/// element type below returns an error rather than an empty rule.
|
||||
///
|
||||
/// The default order is chosen per family to integrate that element's
|
||||
/// stiffness exactly: two Gauss points per direction for tri-/bi-linear
|
||||
/// elements, three for the quadratic ones, and the corresponding
|
||||
/// symmetric rules on simplices.
|
||||
fn quadrature_rule(&self, order: Option<usize>) -> FeaResult<QuadratureRule> {
|
||||
let integration_order = order.unwrap_or(2);
|
||||
// Create quadrature rule based on element type and integration order
|
||||
Ok(quadrature::QuadratureRule::new(
|
||||
vec![],
|
||||
integration_order,
|
||||
self.element_type.spatial_dimension(),
|
||||
))
|
||||
use ElementType::{
|
||||
Hex8, Hex20, Hex27, Line2, Line3, Point, Pyramid5, Pyramid13, Quad4, Quad8, Quad9,
|
||||
Tet4, Tet10, Tri3, Tri6, Wedge6, Wedge15,
|
||||
};
|
||||
|
||||
let order = order.unwrap_or(match self.element_type {
|
||||
Point => 1,
|
||||
Line2 | Tri3 | Tet4 => 1,
|
||||
Line3 | Tri6 | Tet10 | Quad4 | Hex8 | Wedge6 => 2,
|
||||
Quad8 | Quad9 | Hex20 | Hex27 | Wedge15 | Pyramid5 | Pyramid13 => 3,
|
||||
});
|
||||
|
||||
match self.element_type {
|
||||
Point => Ok(QuadratureRule::new(
|
||||
vec![quadrature::QuadraturePoint::new_3d(0.0, 0.0, 0.0, 1.0)],
|
||||
order,
|
||||
0,
|
||||
)),
|
||||
Line2 | Line3 => QuadratureRule::line(order),
|
||||
Tri3 | Tri6 => QuadratureRule::triangle(order),
|
||||
Quad4 | Quad8 | Quad9 => QuadratureRule::quadrilateral(order),
|
||||
Tet4 | Tet10 => QuadratureRule::tetrahedron(order),
|
||||
Hex8 | Hex20 | Hex27 => QuadratureRule::hexahedron(order),
|
||||
|
||||
// A wedge is a triangle extruded along zeta, so its rule is the
|
||||
// tensor product of the triangle rule with a 1-D Gauss rule.
|
||||
Wedge6 | Wedge15 => {
|
||||
let tri = QuadratureRule::triangle(order.min(3))?;
|
||||
let line = QuadratureRule::line(order)?;
|
||||
let mut points = Vec::with_capacity(tri.points.len() * line.points.len());
|
||||
for t in &tri.points {
|
||||
for l in &line.points {
|
||||
points.push(quadrature::QuadraturePoint::new_3d(
|
||||
t.coords.xi(),
|
||||
t.coords.eta(),
|
||||
l.coords.xi(),
|
||||
t.weight * l.weight,
|
||||
));
|
||||
}
|
||||
}
|
||||
Ok(QuadratureRule::new(points, order, 3))
|
||||
}
|
||||
|
||||
// No collapsed-hexahedron rule is implemented for pyramids. An
|
||||
// explicit error beats an empty rule that integrates to zero.
|
||||
Pyramid5 | Pyramid13 => Err(ElementError::UnsupportedQuadratureOrder {
|
||||
element_type: "pyramid".to_string(),
|
||||
requested_order: order,
|
||||
max_order: 0,
|
||||
}
|
||||
.into()),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@@ -653,6 +653,32 @@ impl Quadrilateral9 {
|
||||
pub fn new() -> Self {
|
||||
Self
|
||||
}
|
||||
|
||||
/// Which 1-D Lagrange basis pair `(i, j)` each node uses, with
|
||||
/// `i` indexing `r ∈ {-1, 0, +1}` and `j` indexing `s ∈ {-1, 0, +1}`.
|
||||
///
|
||||
/// The node order is the finite-element convention the rest of this crate
|
||||
/// uses: four corners counter-clockwise, then four mid-edges starting at
|
||||
/// the bottom, then the centre. `Quadrilateral8` and `Quadrilateral4`
|
||||
/// both follow it, and a mesh built for one must be readable by the
|
||||
/// others.
|
||||
///
|
||||
/// Shape functions were previously emitted in raw lexicographic lattice
|
||||
/// order — `(-1,-1), (0,-1), (1,-1), (-1,0), …` — so a mesh written in
|
||||
/// the usual order paired each node with the wrong basis function. At the
|
||||
/// element centre that made the Jacobian exactly singular, because the
|
||||
/// scrambled ordering mapped the element onto a degenerate shape.
|
||||
const NODE_LATTICE: [(usize, usize); 9] = [
|
||||
(0, 0), // 0: (-1, -1)
|
||||
(2, 0), // 1: (+1, -1)
|
||||
(2, 2), // 2: (+1, +1)
|
||||
(0, 2), // 3: (-1, +1)
|
||||
(1, 0), // 4: ( 0, -1)
|
||||
(2, 1), // 5: (+1, 0)
|
||||
(1, 2), // 6: ( 0, +1)
|
||||
(0, 1), // 7: (-1, 0)
|
||||
(1, 1), // 8: ( 0, 0)
|
||||
];
|
||||
}
|
||||
|
||||
impl ShapeFunctions for Quadrilateral9 {
|
||||
@@ -682,13 +708,9 @@ impl ShapeFunctions for Quadrilateral9 {
|
||||
0.5 * s * (s + 1.0), // L₂(s)
|
||||
];
|
||||
|
||||
// 2D tensor product
|
||||
let mut idx = 0;
|
||||
for j in 0..3 {
|
||||
for i in 0..3 {
|
||||
values[idx] = l_r[i] * l_s[j];
|
||||
idx += 1;
|
||||
}
|
||||
// 2D tensor product, emitted in finite-element node order.
|
||||
for (node, &(i, j)) in Self::NODE_LATTICE.iter().enumerate() {
|
||||
values[node] = l_r[i] * l_s[j];
|
||||
}
|
||||
|
||||
Ok(values)
|
||||
@@ -715,14 +737,11 @@ impl ShapeFunctions for Quadrilateral9 {
|
||||
|
||||
let dl_s = [s - 0.5, -2.0 * s, s + 0.5];
|
||||
|
||||
// Compute derivatives
|
||||
let mut idx = 0;
|
||||
for j in 0..3 {
|
||||
for i in 0..3 {
|
||||
derivs[(idx, 0)] = dl_r[i] * l_s[j]; // dN/dr
|
||||
derivs[(idx, 1)] = l_r[i] * dl_s[j]; // dN/ds
|
||||
idx += 1;
|
||||
}
|
||||
// Same node order as `evaluate`, or the Jacobian would pair each
|
||||
// shape function's derivative with a different node's coordinates.
|
||||
for (node, &(i, j)) in Self::NODE_LATTICE.iter().enumerate() {
|
||||
derivs[(node, 0)] = dl_r[i] * l_s[j]; // dN/dr
|
||||
derivs[(node, 1)] = l_r[i] * dl_s[j]; // dN/ds
|
||||
}
|
||||
|
||||
Ok(derivs)
|
||||
|
||||
@@ -593,11 +593,15 @@ impl ShapeFunctions for Hexahedron20 {
|
||||
derivs[(0, 0)] = -0.125 * sm * tm * (-r - s - t - 2.0) - 0.125 * rm * sm * tm;
|
||||
derivs[(1, 0)] = 0.125 * sm * tm * (r - s - t - 2.0) + 0.125 * rp * sm * tm;
|
||||
derivs[(2, 0)] = 0.125 * sp * tm * (r + s - t - 2.0) + 0.125 * rp * sp * tm;
|
||||
derivs[(3, 0)] = -0.125 * sp * tm * (-r + s - t - 2.0) + 0.125 * rm * sp * tm;
|
||||
// N3 = 0.125 rm sp tm (-r + s - t - 2). Both the product rule term
|
||||
// from d(rm)/dr = -1 and the one from d(-r + ...)/dr = -1 are
|
||||
// negative; the second previously carried a `+`.
|
||||
derivs[(3, 0)] = -0.125 * sp * tm * (-r + s - t - 2.0) - 0.125 * rm * sp * tm;
|
||||
derivs[(4, 0)] = -0.125 * sm * tp * (-r - s + t - 2.0) - 0.125 * rm * sm * tp;
|
||||
derivs[(5, 0)] = 0.125 * sm * tp * (r - s + t - 2.0) + 0.125 * rp * sm * tp;
|
||||
derivs[(6, 0)] = 0.125 * sp * tp * (r + s + t - 2.0) + 0.125 * rp * sp * tp;
|
||||
derivs[(7, 0)] = -0.125 * sp * tp * (-r + s + t - 2.0) + 0.125 * rm * sp * tp;
|
||||
// N7 = 0.125 rm sp tp (-r + s + t - 2), same correction as node 3.
|
||||
derivs[(7, 0)] = -0.125 * sp * tp * (-r + s + t - 2.0) - 0.125 * rm * sp * tp;
|
||||
|
||||
// dN/dr for mid-edge nodes
|
||||
derivs[(8, 0)] = -0.5 * r * sm * tm;
|
||||
@@ -616,11 +620,15 @@ impl ShapeFunctions for Hexahedron20 {
|
||||
// dN/ds (similar pattern)
|
||||
// Corner nodes
|
||||
derivs[(0, 1)] = -0.125 * rm * tm * (-r - s - t - 2.0) - 0.125 * rm * sm * tm;
|
||||
derivs[(1, 1)] = -0.125 * rp * tm * (r - s - t - 2.0) + 0.125 * rp * sm * tm;
|
||||
// N1 = 0.125 rp sm tm (r - s - t - 2). d(sm)/ds and d(-s + ...)/ds are
|
||||
// both -1, so both product-rule terms are negative; the second
|
||||
// previously carried a `+`. Same correction as node 5 below, and the
|
||||
// mirror of the one applied to nodes 3 and 7 in dN/dr.
|
||||
derivs[(1, 1)] = -0.125 * rp * tm * (r - s - t - 2.0) - 0.125 * rp * sm * tm;
|
||||
derivs[(2, 1)] = 0.125 * rp * tm * (r + s - t - 2.0) + 0.125 * rp * sp * tm;
|
||||
derivs[(3, 1)] = 0.125 * rm * tm * (-r + s - t - 2.0) + 0.125 * rm * sp * tm;
|
||||
derivs[(4, 1)] = -0.125 * rm * tp * (-r - s + t - 2.0) - 0.125 * rm * sm * tp;
|
||||
derivs[(5, 1)] = -0.125 * rp * tp * (r - s + t - 2.0) + 0.125 * rp * sm * tp;
|
||||
derivs[(5, 1)] = -0.125 * rp * tp * (r - s + t - 2.0) - 0.125 * rp * sm * tp;
|
||||
derivs[(6, 1)] = 0.125 * rp * tp * (r + s + t - 2.0) + 0.125 * rp * sp * tp;
|
||||
derivs[(7, 1)] = 0.125 * rm * tp * (-r + s + t - 2.0) + 0.125 * rm * sp * tp;
|
||||
|
||||
|
||||
@@ -275,11 +275,25 @@ impl ShapeFunctions for Wedge15 {
|
||||
values[6 + i] = tri_funcs[i] * n_plus;
|
||||
}
|
||||
|
||||
// Vertical mid-edges (corners only)
|
||||
// Vertical mid-edges. `4 L_i n₋ n₊` is `L_i (1 - t²)`.
|
||||
values[12] = 4.0 * l1 * n_minus * n_plus; // Vertical edge at corner 1
|
||||
values[13] = 4.0 * l2 * n_minus * n_plus; // Vertical edge at corner 2
|
||||
values[14] = 4.0 * l3 * n_minus * n_plus; // Vertical edge at corner 3
|
||||
|
||||
// Corner correction, without which this element does not form a
|
||||
// partition of unity.
|
||||
//
|
||||
// Introducing a mid-side node on a vertical edge adds `L_i (1 - t²)`
|
||||
// to the sum, so the two corners sharing that edge must each give up
|
||||
// half of it. Omitting the correction left the shape functions
|
||||
// summing to `1 + 4 n₋ n₊`, which is 2 at mid-height — a quadratic
|
||||
// wedge that doubled every field interpolated through it.
|
||||
let bubble = 1.0 - t * t;
|
||||
for (corner, l) in [l1, l2, l3].into_iter().enumerate() {
|
||||
values[corner] -= 0.5 * l * bubble;
|
||||
values[6 + corner] -= 0.5 * l * bubble;
|
||||
}
|
||||
|
||||
Ok(values)
|
||||
}
|
||||
|
||||
@@ -374,6 +388,24 @@ impl ShapeFunctions for Wedge15 {
|
||||
derivs[(13, 2)] = 4.0 * l2 * t_factor; // Vertical edge at corner 2
|
||||
derivs[(14, 2)] = 4.0 * l3 * t_factor; // Vertical edge at corner 3
|
||||
|
||||
// Derivative of the corner correction `-0.5 L_i (1 - t²)` applied in
|
||||
// `evaluate`. The derivatives must carry it too, or they would not be
|
||||
// the derivatives of the functions this element actually uses — and
|
||||
// the sum of dN/dr, dN/ds and dN/dt would not vanish as a partition
|
||||
// of unity requires.
|
||||
let bubble = 1.0 - t * t;
|
||||
// d(l1, l2, l3)/dr and /ds for the area coordinates.
|
||||
let dl_dr = [-1.0, 1.0, 0.0];
|
||||
let dl_ds = [-1.0, 0.0, 1.0];
|
||||
for (corner, l) in [l1, l2, l3].into_iter().enumerate() {
|
||||
for face in [0usize, 6] {
|
||||
derivs[(face + corner, 0)] -= 0.5 * dl_dr[corner] * bubble;
|
||||
derivs[(face + corner, 1)] -= 0.5 * dl_ds[corner] * bubble;
|
||||
// d/dt of `-0.5 L (1 - t²)` is `+L t`.
|
||||
derivs[(face + corner, 2)] += l * t;
|
||||
}
|
||||
}
|
||||
|
||||
Ok(derivs)
|
||||
}
|
||||
|
||||
@@ -808,155 +840,37 @@ impl FiniteElement for Pyramid13 {
|
||||
}
|
||||
|
||||
impl ShapeFunctions for Pyramid13 {
|
||||
fn evaluate(&self, xi: &[f64]) -> FeaResult<DVector<f64>> {
|
||||
if xi.len() < 3 {
|
||||
return Err(crate::error::FeaError::InvalidInput(
|
||||
"Parametric coordinates must have at least 3 components".to_string(),
|
||||
));
|
||||
/// Not implemented.
|
||||
///
|
||||
/// The previous implementation was not a quadratic pyramid basis. Its
|
||||
/// thirteen shape functions summed to 4 at the element centre rather than
|
||||
/// 1, so any field interpolated through it was scaled by a factor that
|
||||
/// varied over the element, and `derivatives` allocated a 13x3 matrix
|
||||
/// then wrote rows 13, 14 and 15 -- it had been copied from a sixteen-node
|
||||
/// layout -- so it panicked with an out-of-bounds index before the wrong
|
||||
/// values could be used.
|
||||
///
|
||||
/// A correct 13-node pyramid basis is rational rather than polynomial,
|
||||
/// and this crate has no pyramid quadrature rule to integrate it with
|
||||
/// (see `StandardFiniteElement::quadrature_rule`), so implementing the
|
||||
/// basis alone would not make the element usable. Both gaps are recorded
|
||||
/// in `omni-cortex/docs/solver_status.md`. Returning an error names the
|
||||
/// limitation at the point of use; `Pyramid5` is unaffected and works.
|
||||
fn evaluate(&self, _xi: &[f64]) -> FeaResult<DVector<f64>> {
|
||||
Err(crate::error::FeaError::InvalidInput(
|
||||
"Pyramid13 shape functions are not implemented; use Pyramid5, or \
|
||||
split the pyramid into tetrahedra"
|
||||
.to_string(),
|
||||
))
|
||||
}
|
||||
|
||||
let r = xi[0];
|
||||
let s = xi[1];
|
||||
let t = xi[2];
|
||||
|
||||
// Pyramid13: quadratic base, linear sides
|
||||
// Nodes 0-3: base corners
|
||||
// Nodes 4-7: base mid-edges
|
||||
// Node 8: base center
|
||||
// Nodes 9-12: vertical mid-edges
|
||||
// Node 13: apex (not used in 13-node version, would be node 12)
|
||||
|
||||
let mut values = DVector::zeros(13);
|
||||
|
||||
// Special handling for degenerate case
|
||||
let eps = 1e-10;
|
||||
let w = if (1.0 - t).abs() < eps { eps } else { 1.0 - t };
|
||||
|
||||
// Modified coordinates for pyramid
|
||||
let _xi_mod = r * w;
|
||||
let _eta_mod = s * w;
|
||||
|
||||
// Base corner nodes
|
||||
values[0] = 0.25 * (1.0 - r) * (1.0 - s) * (1.0 - t) * (-r - s - 2.0 * t + 1.0);
|
||||
values[1] = 0.25 * (1.0 + r) * (1.0 - s) * (1.0 - t) * (r - s - 2.0 * t + 1.0);
|
||||
values[2] = 0.25 * (1.0 + r) * (1.0 + s) * (1.0 - t) * (r + s - 2.0 * t + 1.0);
|
||||
values[3] = 0.25 * (1.0 - r) * (1.0 + s) * (1.0 - t) * (-r + s - 2.0 * t + 1.0);
|
||||
|
||||
// Base mid-edge nodes
|
||||
values[4] = 0.5 * (1.0 - r * r) * (1.0 - s) * (1.0 - t);
|
||||
values[5] = 0.5 * (1.0 + r) * (1.0 - s * s) * (1.0 - t);
|
||||
values[6] = 0.5 * (1.0 - r * r) * (1.0 + s) * (1.0 - t);
|
||||
values[7] = 0.5 * (1.0 - r) * (1.0 - s * s) * (1.0 - t);
|
||||
|
||||
// Base center node
|
||||
values[8] = (1.0 - r * r) * (1.0 - s * s) * (1.0 - t);
|
||||
|
||||
// Vertical mid-edge nodes
|
||||
values[9] = 0.5 * (1.0 - r) * (1.0 - s) * t * (1.0 - t);
|
||||
values[10] = 0.5 * (1.0 + r) * (1.0 - s) * t * (1.0 - t);
|
||||
values[11] = 0.5 * (1.0 + r) * (1.0 + s) * t * (1.0 - t);
|
||||
values[12] = 0.5 * (1.0 - r) * (1.0 + s) * t * (1.0 - t);
|
||||
|
||||
Ok(values)
|
||||
}
|
||||
|
||||
fn derivatives(&self, xi: &[f64]) -> FeaResult<DMatrix<f64>> {
|
||||
if xi.len() < 3 {
|
||||
return Err(crate::error::FeaError::InvalidInput(
|
||||
"Parametric coordinates must have at least 3 components".to_string(),
|
||||
));
|
||||
}
|
||||
|
||||
let r = xi[0];
|
||||
let s = xi[1];
|
||||
let t = xi[2];
|
||||
|
||||
let mut derivs = DMatrix::zeros(13, 3);
|
||||
|
||||
// Simplified derivatives (full implementation would be more complex)
|
||||
// dN/dr for base corner nodes
|
||||
derivs[(0, 0)] = -0.25 * (1.0 - s) * (1.0 - t) * (-2.0 * r - s - 2.0 * t + 1.0);
|
||||
derivs[(1, 0)] = 0.25 * (1.0 - s) * (1.0 - t) * (2.0 * r - s - 2.0 * t + 1.0);
|
||||
derivs[(2, 0)] = 0.25 * (1.0 + s) * (1.0 - t) * (2.0 * r + s - 2.0 * t + 1.0);
|
||||
derivs[(3, 0)] = -0.25 * (1.0 + s) * (1.0 - t) * (-2.0 * r + s - 2.0 * t + 1.0);
|
||||
|
||||
// dN/dr for base mid-edges
|
||||
derivs[(4, 0)] = -r * (1.0 - s) * (1.0 - t);
|
||||
derivs[(5, 0)] = 0.5 * (1.0 - s * s) * (1.0 - t);
|
||||
derivs[(6, 0)] = -r * (1.0 + s) * (1.0 - t);
|
||||
derivs[(7, 0)] = -0.5 * (1.0 - s * s) * (1.0 - t);
|
||||
|
||||
// dN/ds for base corner nodes
|
||||
derivs[(0, 1)] = -0.25 * (1.0 - r) * (1.0 - t) * (-r - 2.0 * s - 2.0 * t + 1.0);
|
||||
derivs[(1, 1)] = -0.25 * (1.0 + r) * (1.0 - t) * (r - 2.0 * s - 2.0 * t + 1.0);
|
||||
derivs[(2, 1)] = 0.25 * (1.0 + r) * (1.0 - t) * (r + 2.0 * s - 2.0 * t + 1.0);
|
||||
derivs[(3, 1)] = 0.25 * (1.0 - r) * (1.0 - t) * (-r + 2.0 * s - 2.0 * t + 1.0);
|
||||
|
||||
// dN/ds for base mid-edges
|
||||
derivs[(4, 1)] = -0.5 * (1.0 - r * r) * (1.0 - t);
|
||||
derivs[(5, 1)] = -s * (1.0 + r) * (1.0 - t);
|
||||
derivs[(6, 1)] = 0.5 * (1.0 - r * r) * (1.0 - t);
|
||||
derivs[(7, 1)] = -s * (1.0 - r) * (1.0 - t);
|
||||
|
||||
// dN/dt for base corner nodes
|
||||
derivs[(0, 2)] = -0.25 * (1.0 - r) * (1.0 - s) * (-r - s - 2.0 * t + 1.0);
|
||||
derivs[(1, 2)] = -0.25 * (1.0 + r) * (1.0 - s) * (r - s - 2.0 * t + 1.0);
|
||||
derivs[(2, 2)] = -0.25 * (1.0 + r) * (1.0 + s) * (r + s - 2.0 * t + 1.0);
|
||||
derivs[(3, 2)] = -0.25 * (1.0 - r) * (1.0 + s) * (-r + s - 2.0 * t + 1.0);
|
||||
|
||||
// dN/dt for base mid-edges
|
||||
derivs[(4, 2)] = -0.5 * (1.0 - r * r) * (1.0 - s);
|
||||
derivs[(5, 2)] = -0.5 * (1.0 - s * s) * (1.0 + r);
|
||||
derivs[(6, 2)] = -0.5 * (1.0 - r * r) * (1.0 + s);
|
||||
derivs[(7, 2)] = -0.5 * (1.0 - s * s) * (1.0 - r);
|
||||
|
||||
// Top layer nodes (8-15) - similar to base but with (1+t) factor
|
||||
// dN/dr for top corner nodes
|
||||
derivs[(8, 0)] = -0.25 * (1.0 - s) * (1.0 + t) * (-2.0 * r - s + 2.0 * t + 1.0);
|
||||
derivs[(9, 0)] = 0.25 * (1.0 - s) * (1.0 + t) * (2.0 * r - s + 2.0 * t + 1.0);
|
||||
derivs[(10, 0)] = 0.25 * (1.0 + s) * (1.0 + t) * (2.0 * r + s + 2.0 * t + 1.0);
|
||||
derivs[(11, 0)] = -0.25 * (1.0 + s) * (1.0 + t) * (-2.0 * r + s + 2.0 * t + 1.0);
|
||||
|
||||
// dN/dr for top mid-edges
|
||||
derivs[(12, 0)] = -r * (1.0 - s) * (1.0 + t);
|
||||
derivs[(13, 0)] = 0.5 * (1.0 - s * s) * (1.0 + t);
|
||||
derivs[(14, 0)] = -r * (1.0 + s) * (1.0 + t);
|
||||
derivs[(15, 0)] = -0.5 * (1.0 - s * s) * (1.0 + t);
|
||||
|
||||
// dN/ds for top corner nodes
|
||||
derivs[(8, 1)] = -0.25 * (1.0 - r) * (1.0 + t) * (-r - 2.0 * s + 2.0 * t + 1.0);
|
||||
derivs[(9, 1)] = -0.25 * (1.0 + r) * (1.0 + t) * (r - 2.0 * s + 2.0 * t + 1.0);
|
||||
derivs[(10, 1)] = 0.25 * (1.0 + r) * (1.0 + t) * (r + 2.0 * s + 2.0 * t + 1.0);
|
||||
derivs[(11, 1)] = 0.25 * (1.0 - r) * (1.0 + t) * (-r + 2.0 * s + 2.0 * t + 1.0);
|
||||
|
||||
// dN/ds for top mid-edges
|
||||
derivs[(12, 1)] = -0.5 * (1.0 - r * r) * (1.0 + t);
|
||||
derivs[(13, 1)] = -s * (1.0 + r) * (1.0 + t);
|
||||
derivs[(14, 1)] = 0.5 * (1.0 - r * r) * (1.0 + t);
|
||||
derivs[(15, 1)] = -s * (1.0 - r) * (1.0 + t);
|
||||
|
||||
// dN/dt for top corner nodes
|
||||
derivs[(8, 2)] = 0.25 * (1.0 - r) * (1.0 - s) * (-r - s + 2.0 * t + 1.0);
|
||||
derivs[(9, 2)] = 0.25 * (1.0 + r) * (1.0 - s) * (r - s + 2.0 * t + 1.0);
|
||||
derivs[(10, 2)] = 0.25 * (1.0 + r) * (1.0 + s) * (r + s + 2.0 * t + 1.0);
|
||||
derivs[(11, 2)] = 0.25 * (1.0 - r) * (1.0 + s) * (-r + s + 2.0 * t + 1.0);
|
||||
|
||||
// dN/dt for top mid-edges
|
||||
derivs[(12, 2)] = 0.5 * (1.0 - r * r) * (1.0 - s);
|
||||
derivs[(13, 2)] = 0.5 * (1.0 - s * s) * (1.0 + r);
|
||||
derivs[(14, 2)] = 0.5 * (1.0 - r * r) * (1.0 + s);
|
||||
derivs[(15, 2)] = 0.5 * (1.0 - s * s) * (1.0 - r);
|
||||
|
||||
// Vertical mid-edges (16-19)
|
||||
for i in 16..20 {
|
||||
let _base_idx = i - 16; // 0,1,2,3
|
||||
// These are the derivatives for vertical edge shape functions
|
||||
derivs[(i, 0)] = 0.0; // No r-derivative for vertical edges
|
||||
derivs[(i, 1)] = 0.0; // No s-derivative for vertical edges
|
||||
derivs[(i, 2)] = 1.0; // Unit t-derivative for vertical edges
|
||||
}
|
||||
|
||||
Ok(derivs)
|
||||
/// Not implemented. See [`Pyramid13::evaluate`].
|
||||
fn derivatives(&self, _xi: &[f64]) -> FeaResult<DMatrix<f64>> {
|
||||
Err(crate::error::FeaError::InvalidInput(
|
||||
"Pyramid13 shape function derivatives are not implemented; use \
|
||||
Pyramid5, or split the pyramid into tetrahedra"
|
||||
.to_string(),
|
||||
))
|
||||
}
|
||||
|
||||
fn num_nodes(&self) -> usize {
|
||||
@@ -1024,12 +938,17 @@ mod tests {
|
||||
|
||||
#[test]
|
||||
fn test_pyramid13() {
|
||||
let elem = Pyramid13;
|
||||
let xi = vec![0.0, 0.0, 0.0];
|
||||
let shape = elem.evaluate(&xi).unwrap();
|
||||
|
||||
assert_eq!(shape.len(), 13);
|
||||
// Note: Sum may not be exactly 1.0 due to the special pyramid formulation
|
||||
assert!(shape.sum() > 0.0);
|
||||
// Pyramid13 is not implemented and reports so. See the doc comment on
|
||||
// its `evaluate`.
|
||||
//
|
||||
// This test previously accepted any positive sum, explaining that the
|
||||
// "special pyramid formulation" need not sum to 1. No formulation is
|
||||
// exempt: shape functions that do not sum to unity cannot reproduce a
|
||||
// constant field, so the element could not represent a rigid-body
|
||||
// translation. The actual sum was 4 at the element centre.
|
||||
let error = Pyramid13
|
||||
.evaluate(&[0.0, 0.0, 0.0])
|
||||
.expect_err("Pyramid13 must report that it is unimplemented");
|
||||
assert!(error.to_string().to_lowercase().contains("not implemented"));
|
||||
}
|
||||
}
|
||||
|
||||
@@ -287,10 +287,14 @@ pub struct MemoryInfo {
|
||||
|
||||
impl std::fmt::Display for MemoryInfo {
|
||||
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
|
||||
// Binary units. Driver APIs report memory in bytes and the sizes are
|
||||
// powers of two, so 1 GiB = 1024^3; dividing by 1e9 instead labelled
|
||||
// an 8 GiB device "8.59 GB".
|
||||
const GIB: f64 = (1024 * 1024 * 1024) as f64;
|
||||
writeln!(f, "GPU Memory Information:")?;
|
||||
writeln!(f, " Total: {:.2} GB", self.total_memory as f64 / 1e9)?;
|
||||
writeln!(f, " Free: {:.2} GB", self.free_memory as f64 / 1e9)?;
|
||||
writeln!(f, " Used: {:.2} GB", self.used_memory as f64 / 1e9)?;
|
||||
writeln!(f, " Total: {:.2} GiB", self.total_memory as f64 / GIB)?;
|
||||
writeln!(f, " Free: {:.2} GiB", self.free_memory as f64 / GIB)?;
|
||||
writeln!(f, " Used: {:.2} GiB", self.used_memory as f64 / GIB)?;
|
||||
writeln!(
|
||||
f,
|
||||
" Usage: {:.1}%",
|
||||
@@ -585,11 +589,14 @@ mod tests {
|
||||
used_memory: 2 * 1024 * 1024 * 1024, // 2 GB
|
||||
};
|
||||
|
||||
// The fixture builds binary gibibytes, so the display must report
|
||||
// them as such. This previously expected "8.00 GB" from a formatter
|
||||
// that divided by 1e9, which renders 8 GiB as 8.59.
|
||||
let display_str = format!("{}", info);
|
||||
assert!(display_str.contains("Total: 8.00 GB"));
|
||||
assert!(display_str.contains("Free: 6.00 GB"));
|
||||
assert!(display_str.contains("Used: 2.00 GB"));
|
||||
assert!(display_str.contains("Usage: 25.0%"));
|
||||
assert!(display_str.contains("Total: 8.00 GiB"), "{display_str}");
|
||||
assert!(display_str.contains("Free: 6.00 GiB"), "{display_str}");
|
||||
assert!(display_str.contains("Used: 2.00 GiB"), "{display_str}");
|
||||
assert!(display_str.contains("Usage: 25.0%"), "{display_str}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
|
||||
@@ -14,29 +14,68 @@
|
||||
//! - **Efficient Assembly**: Sparse matrix assembly with GPU optimization
|
||||
//! - **Robust Solvers**: Direct and iterative solvers with GPU acceleration
|
||||
//! - **Mesh Management**: Advanced mesh operations including refinement and partitioning
|
||||
//! - **Production Ready**: No mocks, stubs, or TODOs - complete implementation
|
||||
//!
|
||||
//! ## Maturity
|
||||
//!
|
||||
//! Validated: element stiffness and mass matrices, shape functions across the
|
||||
//! element library, global assembly, constraint handling, the generalized
|
||||
//! eigensolver, and modal analysis end to end against closed-form bar
|
||||
//! frequencies. See `tests/element_matrices_physical.rs`,
|
||||
//! `tests/shape_function_invariants.rs`, `tests/eigenvalue_closed_form.rs`
|
||||
//! and `tests/modal_closed_form.rs`.
|
||||
//!
|
||||
//! Not yet: `DynamicAnalysis` and `NonlinearAnalysis` return zeros;
|
||||
//! `Pyramid13` and pyramid quadrature are unimplemented and report so. Much
|
||||
//! of the crate's test suite is still disabled behind `#[cfg(disabled)]`, so
|
||||
//! a green run does not imply coverage of what those modules cover.
|
||||
//!
|
||||
//! ## Quick Start
|
||||
//!
|
||||
//! ```rust,no_run
|
||||
//! use rtx_fea::{
|
||||
//! mesh::{Mesh, geometry::Rectangle},
|
||||
//! materials::LinearElastic,
|
||||
//! solvers::DirectSolver,
|
||||
//! analysis::StaticAnalysis,
|
||||
//! Natural frequencies of a fixed-free bar:
|
||||
//!
|
||||
//! ```rust
|
||||
//! use rtx_fea::analysis::{Analysis, AnalysisConfig, AnalysisData, ModalAnalysis};
|
||||
//! use rtx_fea::assembly::DofComponent;
|
||||
//! use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, DirichletBC};
|
||||
//! use rtx_fea::materials::{LinearElastic, MaterialDatabase};
|
||||
//! use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node};
|
||||
//!
|
||||
//! // A 4 x 1 grid of quadrilaterals spanning a 1.0 x 0.05 strip.
|
||||
//! let mut mesh = Mesh::new(2)?;
|
||||
//! let mut columns = Vec::new();
|
||||
//! for i in 0..=4 {
|
||||
//! let x = f64::from(i) * 0.25;
|
||||
//! columns.push([
|
||||
//! mesh.add_node(Node::new_2d(x, 0.0)),
|
||||
//! mesh.add_node(Node::new_2d(x, 0.05)),
|
||||
//! ]);
|
||||
//! }
|
||||
//! for i in 0..4 {
|
||||
//! let nodes = vec![columns[i][0], columns[i + 1][0], columns[i + 1][1], columns[i][1]];
|
||||
//! mesh.add_element(Element::new(ElementType::Quad4, nodes, MaterialId(0))?)?;
|
||||
//! }
|
||||
//!
|
||||
//! let mut materials = MaterialDatabase::new();
|
||||
//! materials.add_material(MaterialId(0), LinearElastic::new(200e9, 0.3).with_density(8000.0), None);
|
||||
//!
|
||||
//! // Clamp the left edge. Without constraints the stiffness matrix is
|
||||
//! // singular and the analysis reports that rather than returning noise.
|
||||
//! let mut boundary_conditions = BoundaryConditionSet::new();
|
||||
//! boundary_conditions.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
|
||||
//! columns[0].to_vec(),
|
||||
//! vec![DofComponent::DisplacementX, DofComponent::DisplacementY],
|
||||
//! 0.0,
|
||||
//! )));
|
||||
//!
|
||||
//! let mut analysis = ModalAnalysis::new(mesh, materials, 3, AnalysisConfig::default())
|
||||
//! .with_boundary_conditions(boundary_conditions);
|
||||
//! let results = analysis.run()?;
|
||||
//!
|
||||
//! let AnalysisData::Vector(frequencies) = &results.additional_data["frequencies"] else {
|
||||
//! panic!("frequencies should be a vector")
|
||||
//! };
|
||||
//!
|
||||
//! // Create a mesh
|
||||
//! let rect = Rectangle::new(1.0, 1.0);
|
||||
//! let mesh = rect.generate_quad_mesh(10, 10, rtx_fea::mesh::MaterialId(0))?;
|
||||
//!
|
||||
//! // Set up analysis
|
||||
//! let material = LinearElastic::new(200e9, 0.3); // Steel
|
||||
//! let mut analysis = StaticAnalysis::new(mesh);
|
||||
//! analysis.add_material(material);
|
||||
//!
|
||||
//! // Apply boundary conditions and solve
|
||||
//! // ... (see examples for complete workflow)
|
||||
//! assert_eq!(frequencies.len(), 3);
|
||||
//! assert!(frequencies.iter().all(|f| f.is_finite() && *f > 0.0));
|
||||
//! # Ok::<(), Box<dyn std::error::Error>>(())
|
||||
//! ```
|
||||
//!
|
||||
|
||||
@@ -24,6 +24,7 @@ use cudarc::driver::safe::CudaContext;
|
||||
use nalgebra::{DMatrix, Vector3, Vector6};
|
||||
use serde::{Deserialize, Serialize};
|
||||
use std::collections::HashMap;
|
||||
use std::sync::Arc;
|
||||
|
||||
pub use composite::*;
|
||||
pub use damage::*;
|
||||
@@ -333,7 +334,7 @@ pub trait Material: Send + Sync {
|
||||
/// Material database for managing multiple materials.
|
||||
pub struct MaterialDatabase {
|
||||
/// Materials storage
|
||||
materials: HashMap<crate::mesh::MaterialId, Box<dyn Material>>,
|
||||
materials: HashMap<crate::mesh::MaterialId, Arc<dyn Material>>,
|
||||
/// Material names
|
||||
names: HashMap<crate::mesh::MaterialId, String>,
|
||||
}
|
||||
@@ -364,7 +365,7 @@ impl MaterialDatabase {
|
||||
name: Option<String>,
|
||||
) {
|
||||
let name = name.unwrap_or_else(|| format!("Material_{}", id.as_usize()));
|
||||
self.materials.insert(id, Box::new(material));
|
||||
self.materials.insert(id, Arc::new(material));
|
||||
self.names.insert(id, name);
|
||||
}
|
||||
|
||||
@@ -403,16 +404,23 @@ impl MaterialDatabase {
|
||||
}
|
||||
|
||||
impl Clone for MaterialDatabase {
|
||||
/// Clone the database, sharing the materials themselves.
|
||||
///
|
||||
/// Materials are immutable once registered, so the `Arc` handles can be
|
||||
/// shared rather than deep-copied and no `clone_box` on the trait object
|
||||
/// is needed.
|
||||
///
|
||||
/// This previously cloned **names only** and silently dropped every
|
||||
/// material, on the grounds that `Box<dyn Material>` cannot be cloned.
|
||||
/// Because `GlobalAssembler` is constructed with `materials.clone()`,
|
||||
/// that meant every assembler ever built received an empty database and
|
||||
/// every analysis failed with `MaterialNotFound` — for a mesh that was
|
||||
/// correctly specified.
|
||||
fn clone(&self) -> Self {
|
||||
// For now, create an empty database with same names
|
||||
// Material cloning would require dyn trait object cloning support
|
||||
let mut new_db = Self::new();
|
||||
for (id, name) in &self.names {
|
||||
new_db.names.insert(*id, name.clone());
|
||||
Self {
|
||||
materials: self.materials.clone(),
|
||||
names: self.names.clone(),
|
||||
}
|
||||
// Note: materials are not cloned, this is a shallow clone for names only
|
||||
// Full implementation would require Material trait to have clone_box method
|
||||
new_db
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@@ -102,6 +102,36 @@ impl EigenvalueSolver {
|
||||
let x = solve_lower(&l, &shifted)?;
|
||||
let b = solve_lower(&l, &x.transpose())?.transpose();
|
||||
|
||||
// Reject a rank-deficient problem explicitly.
|
||||
//
|
||||
// `try_inverse` does not fail on a matrix that is singular only to
|
||||
// working precision, and an unconstrained structure is exactly that
|
||||
// case: its rigid-body modes are zero eigenvalues that round to
|
||||
// ~1e-16 of the largest. Inverting anyway succeeds and returns
|
||||
// enormous, meaningless entries, which emerge as a handful of
|
||||
// near-zero frequencies that look like real low-frequency modes.
|
||||
// Better to say what is wrong and name the remedy.
|
||||
//
|
||||
// Only meaningful without a shift: `K - σM` is deliberately
|
||||
// indefinite when a shift is used to bracket interior modes.
|
||||
if self.shift.is_none() {
|
||||
let spectrum = b.clone().symmetric_eigenvalues();
|
||||
let largest = spectrum.iter().fold(0.0_f64, |acc, v| acc.max(v.abs()));
|
||||
let smallest = spectrum
|
||||
.iter()
|
||||
.fold(f64::INFINITY, |acc: f64, v| acc.min(v.abs()));
|
||||
|
||||
if largest <= 0.0 || smallest / largest < 1e-12 {
|
||||
return Err(SolverError::SolveError {
|
||||
reason: "stiffness matrix is singular — an unconstrained structure has \
|
||||
rigid-body modes at zero frequency. Constrain it, or pass an \
|
||||
explicit shift via `with_shift` to compute the modes around one."
|
||||
.to_string(),
|
||||
}
|
||||
.into());
|
||||
}
|
||||
}
|
||||
|
||||
// A = B⁻¹. Symmetrized to shed the asymmetry that the inverse of a
|
||||
// numerically-symmetric matrix picks up at rounding level; Lanczos is
|
||||
// sensitive to it.
|
||||
|
||||
@@ -12,7 +12,20 @@ mod tests {
|
||||
fn test_math_utils_constants() {
|
||||
assert!(MathUtils::EPSILON > 0.0);
|
||||
assert!(MathUtils::SMALL > 0.0);
|
||||
assert!(MathUtils::SMALL < MathUtils::EPSILON * 1000.0);
|
||||
|
||||
// SMALL is the practical "is this zero" threshold and must sit well
|
||||
// above machine epsilon, or comparisons against it would be decided
|
||||
// by rounding noise.
|
||||
//
|
||||
// This previously asserted `SMALL < EPSILON * 1000`, which is false
|
||||
// and inverts the relationship: EPSILON is f64::EPSILON (2.2e-16), so
|
||||
// that bound is 2.2e-13, while SMALL is 1e-12. A threshold below a
|
||||
// thousand machine epsilons would defeat its own purpose.
|
||||
assert!(MathUtils::SMALL > MathUtils::EPSILON);
|
||||
|
||||
// ...and far below any length, stress or stiffness a model uses, so
|
||||
// it never swallows a physically meaningful quantity.
|
||||
assert!(MathUtils::SMALL < 1e-6);
|
||||
}
|
||||
|
||||
#[test]
|
||||
@@ -198,16 +211,31 @@ mod tests {
|
||||
|
||||
#[test]
|
||||
fn test_von_mises_stress_biaxial() {
|
||||
// Equal biaxial stress
|
||||
// Equal biaxial stress: sigma_x = sigma_y = 100, sigma_z = 0.
|
||||
//
|
||||
// sigma_vm = sqrt(1/2 [(sx-sy)^2 + (sy-sz)^2 + (sz-sx)^2])
|
||||
// = sqrt(1/2 [0 + 100^2 + 100^2]) = 100.
|
||||
//
|
||||
// This previously expected 0, commented "no deviatoric stress". That
|
||||
// is only true of a *hydrostatic* state, where all three principal
|
||||
// stresses are equal. An equal biaxial state with a free third
|
||||
// direction is strongly deviatoric -- and expecting zero would mean a
|
||||
// biaxially loaded sheet could never yield.
|
||||
let stress = DVector::from_vec(vec![100.0, 100.0, 0.0, 0.0, 0.0, 0.0]);
|
||||
let vm_stress = VoigtUtils::von_mises_stress(&stress);
|
||||
assert!(MathUtils::approx_eq(vm_stress, 0.0, 1e-10)); // No deviatoric stress
|
||||
assert!(MathUtils::approx_eq(vm_stress, 100.0, 1e-10));
|
||||
|
||||
// Unequal biaxial stress
|
||||
// Genuine hydrostatic stress, which does have zero von Mises stress.
|
||||
// This is the case the comment above was reaching for.
|
||||
let stress = DVector::from_vec(vec![100.0, 100.0, 100.0, 0.0, 0.0, 0.0]);
|
||||
let vm_stress = VoigtUtils::von_mises_stress(&stress);
|
||||
assert!(MathUtils::approx_eq(vm_stress, 0.0, 1e-10));
|
||||
|
||||
// Unequal biaxial stress: sqrt(1/2 [50^2 + 50^2 + 100^2]) = sqrt(7500).
|
||||
// Not |100 - 50|; the von Mises stress is not a principal difference.
|
||||
let stress = DVector::from_vec(vec![100.0, 50.0, 0.0, 0.0, 0.0, 0.0]);
|
||||
let vm_stress = VoigtUtils::von_mises_stress(&stress);
|
||||
let expected = 50.0; // |100 - 50| = 50
|
||||
assert!(MathUtils::approx_eq(vm_stress, expected, 1e-10));
|
||||
assert!(MathUtils::approx_eq(vm_stress, 7500.0_f64.sqrt(), 1e-10));
|
||||
}
|
||||
|
||||
#[test]
|
||||
@@ -281,13 +309,43 @@ mod tests {
|
||||
fn test_integrate_1d_trigonometric() {
|
||||
use std::f64::consts::PI;
|
||||
|
||||
// Test ∫ sin(x) dx from 0 to π = 2
|
||||
let result = IntegrationUtils::integrate_1d(|x| x.sin(), 0.0, PI, 3).unwrap();
|
||||
assert!(MathUtils::approx_eq(result, 2.0, 1e-10));
|
||||
// Gauss-Legendre with n points is exact for polynomials of degree
|
||||
// 2n-1, and sin is not a polynomial. Three points therefore give an
|
||||
// approximation, not an exact answer -- the error here is ~1e-3.
|
||||
//
|
||||
// This previously demanded 1e-10 from a 3-point rule, which no
|
||||
// correct implementation could satisfy. Rather than relax the bound
|
||||
// to whatever happens to pass, assert what actually pins the
|
||||
// quadrature: that the error *shrinks* as points are added. A wrong
|
||||
// rule can land near the right answer once; it will not converge.
|
||||
let exact = 2.0;
|
||||
let errors: Vec<f64> = (1..=3)
|
||||
.map(|n| {
|
||||
let result = IntegrationUtils::integrate_1d(|x| x.sin(), 0.0, PI, n).unwrap();
|
||||
(result - exact).abs()
|
||||
})
|
||||
.collect();
|
||||
|
||||
// Test ∫ cos(x) dx from 0 to π/2 = 1
|
||||
for pair in errors.windows(2) {
|
||||
assert!(
|
||||
pair[1] < pair[0],
|
||||
"adding a Gauss point did not reduce the error: {errors:?}"
|
||||
);
|
||||
}
|
||||
assert!(
|
||||
errors[2] < 2e-3,
|
||||
"3-point Gauss on sin over [0, pi] should be accurate to ~1e-3, got {}",
|
||||
errors[2]
|
||||
);
|
||||
|
||||
// Same for cos over a quarter period.
|
||||
let result = IntegrationUtils::integrate_1d(|x| x.cos(), 0.0, PI / 2.0, 3).unwrap();
|
||||
assert!(MathUtils::approx_eq(result, 1.0, 1e-10));
|
||||
assert!(MathUtils::approx_eq(result, 1.0, 1e-3));
|
||||
|
||||
// A cubic *is* within the exactness of a 2-point rule, so there the
|
||||
// tight tolerance is the right expectation.
|
||||
let result = IntegrationUtils::integrate_1d(|x| x * x * x, 0.0, 2.0, 2).unwrap();
|
||||
assert!(MathUtils::approx_eq(result, 4.0, 1e-12));
|
||||
}
|
||||
|
||||
#[test]
|
||||
|
||||
@@ -0,0 +1,290 @@
|
||||
//! Physical invariants of the element stiffness and mass matrices.
|
||||
//!
|
||||
//! These are the cheapest decisive checks that element matrix computation is
|
||||
//! real. None of them needs a reference table: each is a property the exact
|
||||
//! matrices must satisfy for any correct implementation, so a failure
|
||||
//! localises to the quadrature, the Jacobian, the shape function derivatives
|
||||
//! or the constitutive matrix rather than to an accuracy budget.
|
||||
//!
|
||||
//! They exist because `StandardFiniteElement::compute_element_matrices`
|
||||
//! returned `DMatrix::zeros(..)` for every element in the mesh, which made
|
||||
//! every global matrix `GlobalAssembler` produced zero, and every analysis
|
||||
//! built on it a solve of a null system. A zero matrix passes a symmetry
|
||||
//! check and a "does not crash" check, so the invariants below are chosen to
|
||||
//! be ones a zero matrix fails.
|
||||
|
||||
use approx::assert_relative_eq;
|
||||
use nalgebra::{DVector, Vector3};
|
||||
use rtx_fea::elements::StandardFiniteElement;
|
||||
use rtx_fea::materials::LinearElastic;
|
||||
use rtx_fea::mesh::ElementType;
|
||||
|
||||
const E: f64 = 210e9;
|
||||
const NU: f64 = 0.3;
|
||||
const RHO: f64 = 7850.0;
|
||||
|
||||
fn steel() -> LinearElastic {
|
||||
LinearElastic::new(E, NU).with_density(RHO)
|
||||
}
|
||||
|
||||
/// Unit square Quad4, counter-clockwise from the origin. Area = 1.
|
||||
fn unit_square() -> StandardFiniteElement {
|
||||
StandardFiniteElement::new(
|
||||
ElementType::Quad4,
|
||||
vec![
|
||||
Vector3::new(0.0, 0.0, 0.0),
|
||||
Vector3::new(1.0, 0.0, 0.0),
|
||||
Vector3::new(1.0, 1.0, 0.0),
|
||||
Vector3::new(0.0, 1.0, 0.0),
|
||||
],
|
||||
)
|
||||
}
|
||||
|
||||
/// A 2 x 3 rectangle, to catch anything that only works when `det J == 1`.
|
||||
fn rectangle_2x3() -> StandardFiniteElement {
|
||||
StandardFiniteElement::new(
|
||||
ElementType::Quad4,
|
||||
vec![
|
||||
Vector3::new(0.0, 0.0, 0.0),
|
||||
Vector3::new(2.0, 0.0, 0.0),
|
||||
Vector3::new(2.0, 3.0, 0.0),
|
||||
Vector3::new(0.0, 3.0, 0.0),
|
||||
],
|
||||
)
|
||||
}
|
||||
|
||||
/// The stub this suite exists to catch returned zeros everywhere.
|
||||
#[test]
|
||||
fn element_matrices_are_not_zero() {
|
||||
let m = unit_square()
|
||||
.compute_element_matrices(&steel(), 0.0)
|
||||
.unwrap();
|
||||
|
||||
assert!(
|
||||
m.stiffness_matrix.amax() > 0.0,
|
||||
"stiffness matrix is entirely zero — element matrix computation is a stub"
|
||||
);
|
||||
let mass = m
|
||||
.mass_matrix
|
||||
.as_ref()
|
||||
.expect("a mass matrix is required for modal and dynamic analysis");
|
||||
assert!(
|
||||
mass.amax() > 0.0,
|
||||
"mass matrix is entirely zero — element matrix computation is a stub"
|
||||
);
|
||||
}
|
||||
|
||||
/// `K` must be symmetric: it is `∫ Bᵀ D B dV` with `D` symmetric.
|
||||
#[test]
|
||||
fn stiffness_is_symmetric() {
|
||||
let m = rectangle_2x3()
|
||||
.compute_element_matrices(&steel(), 0.0)
|
||||
.unwrap();
|
||||
let k = &m.stiffness_matrix;
|
||||
|
||||
for i in 0..k.nrows() {
|
||||
for j in 0..k.ncols() {
|
||||
assert_relative_eq!(k[(i, j)], k[(j, i)], epsilon = 1e-6 * k.amax());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// A rigid-body translation stores no strain energy, so `K t = 0`.
|
||||
///
|
||||
/// This is the single most informative check on the strain-displacement
|
||||
/// matrix `B`: it fails for a sign error, a mis-transformed derivative, or a
|
||||
/// wrong DOF ordering, none of which a symmetry check detects. It is the
|
||||
/// element-level half of the patch test.
|
||||
#[test]
|
||||
fn rigid_translation_produces_no_internal_force() {
|
||||
let element = rectangle_2x3();
|
||||
let m = element.compute_element_matrices(&steel(), 0.0).unwrap();
|
||||
let k = &m.stiffness_matrix;
|
||||
|
||||
let num_nodes = 4;
|
||||
let dim = 2;
|
||||
assert_eq!(k.nrows(), num_nodes * dim);
|
||||
|
||||
// One translation per spatial direction.
|
||||
for d in 0..dim {
|
||||
let mut t = DVector::zeros(num_nodes * dim);
|
||||
for node in 0..num_nodes {
|
||||
t[node * dim + d] = 1.0;
|
||||
}
|
||||
|
||||
let f = k * &t;
|
||||
assert!(
|
||||
f.amax() < 1e-6 * k.amax(),
|
||||
"translating the element in direction {d} produced internal force {:.3e} \
|
||||
against a stiffness scale of {:.3e}",
|
||||
f.amax(),
|
||||
k.amax()
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// An unconstrained plane element has exactly three rigid-body modes — two
|
||||
/// translations and one rotation — so `K` has a three-dimensional null space.
|
||||
///
|
||||
/// Fewer means the element is spuriously stiff; more means it is rank
|
||||
/// deficient and admits a zero-energy deformation (hourglassing), which shows
|
||||
/// up in a real analysis as a mode shape that is pure noise.
|
||||
#[test]
|
||||
fn stiffness_has_exactly_three_rigid_body_modes_in_2d() {
|
||||
let m = rectangle_2x3()
|
||||
.compute_element_matrices(&steel(), 0.0)
|
||||
.unwrap();
|
||||
let k = m.stiffness_matrix;
|
||||
|
||||
let eigenvalues = k.clone().symmetric_eigenvalues();
|
||||
let scale = k.amax();
|
||||
let num_zero = eigenvalues
|
||||
.iter()
|
||||
.filter(|&&e| e.abs() < 1e-9 * scale)
|
||||
.count();
|
||||
|
||||
assert_eq!(
|
||||
num_zero,
|
||||
3,
|
||||
"expected 3 rigid-body modes (2 translations + 1 rotation), found {num_zero}; \
|
||||
eigenvalues (scaled): {:?}",
|
||||
eigenvalues.iter().map(|e| e / scale).collect::<Vec<_>>()
|
||||
);
|
||||
}
|
||||
|
||||
/// The consistent mass matrix integrates to the element's total mass.
|
||||
///
|
||||
/// `Σᵢⱼ Mᵢⱼ = ∫ ρ (Σᵢ Nᵢ)(Σⱼ Nⱼ) dV = ∫ ρ dV = ρV`, using the partition of
|
||||
/// unity. This is decisive for the quadrature rule and the Jacobian
|
||||
/// determinant together: get either wrong and the total mass is wrong by
|
||||
/// exactly the factor of the error.
|
||||
#[test]
|
||||
fn consistent_mass_integrates_to_rho_times_volume() {
|
||||
for (element, area) in [(unit_square(), 1.0), (rectangle_2x3(), 6.0)] {
|
||||
let m = element.compute_element_matrices(&steel(), 0.0).unwrap();
|
||||
let mass = m.mass_matrix.expect("mass matrix required");
|
||||
|
||||
// Vector-valued mass: each spatial direction carries the full mass,
|
||||
// so the whole matrix sums to `dim * rho * V`.
|
||||
let dim = 2;
|
||||
let total: f64 = mass.iter().sum();
|
||||
assert_relative_eq!(total, dim as f64 * RHO * area, max_relative = 1e-9);
|
||||
}
|
||||
}
|
||||
|
||||
/// The mass matrix must be positive definite.
|
||||
///
|
||||
/// It is the metric in `K φ = λ M φ`; if it is not positive definite the
|
||||
/// eigenproblem has no real spectrum and the Cholesky reduction in the
|
||||
/// eigensolver fails outright.
|
||||
#[test]
|
||||
fn mass_matrix_is_positive_definite() {
|
||||
let m = rectangle_2x3()
|
||||
.compute_element_matrices(&steel(), 0.0)
|
||||
.unwrap();
|
||||
let mass = m.mass_matrix.expect("mass matrix required");
|
||||
|
||||
let eigenvalues = mass.clone().symmetric_eigenvalues();
|
||||
let min = eigenvalues.iter().cloned().fold(f64::INFINITY, f64::min);
|
||||
assert!(
|
||||
min > 0.0,
|
||||
"mass matrix has a non-positive eigenvalue {min:.3e}; \
|
||||
eigenvalues: {eigenvalues:?}"
|
||||
);
|
||||
}
|
||||
|
||||
/// Mass scales linearly with density, stiffness with elastic modulus, and
|
||||
/// neither picks up the other's property.
|
||||
///
|
||||
/// A single hard-coded default leaking into the computation would break this,
|
||||
/// and would otherwise be invisible in a suite that only ever uses one
|
||||
/// material.
|
||||
#[test]
|
||||
fn matrices_scale_with_the_material_they_are_given() {
|
||||
let base = steel();
|
||||
let stiffer = LinearElastic::new(2.0 * E, NU).with_density(RHO);
|
||||
let denser = LinearElastic::new(E, NU).with_density(3.0 * RHO);
|
||||
|
||||
let element = rectangle_2x3();
|
||||
let m0 = element.compute_element_matrices(&base, 0.0).unwrap();
|
||||
let m_stiff = element.compute_element_matrices(&stiffer, 0.0).unwrap();
|
||||
let m_dense = element.compute_element_matrices(&denser, 0.0).unwrap();
|
||||
|
||||
assert_relative_eq!(
|
||||
m_stiff.stiffness_matrix.amax(),
|
||||
2.0 * m0.stiffness_matrix.amax(),
|
||||
max_relative = 1e-12
|
||||
);
|
||||
assert_relative_eq!(
|
||||
m_dense.mass_matrix.as_ref().unwrap().amax(),
|
||||
3.0 * m0.mass_matrix.as_ref().unwrap().amax(),
|
||||
max_relative = 1e-12
|
||||
);
|
||||
|
||||
// Doubling E must not change the mass; tripling ρ must not change K.
|
||||
assert_relative_eq!(
|
||||
m_stiff.mass_matrix.as_ref().unwrap().amax(),
|
||||
m0.mass_matrix.as_ref().unwrap().amax(),
|
||||
max_relative = 1e-12
|
||||
);
|
||||
assert_relative_eq!(
|
||||
m_dense.stiffness_matrix.amax(),
|
||||
m0.stiffness_matrix.amax(),
|
||||
max_relative = 1e-12
|
||||
);
|
||||
}
|
||||
|
||||
/// The same invariants must hold in 3D, where the null space is six
|
||||
/// dimensional — three translations and three rotations.
|
||||
#[test]
|
||||
fn hexahedron_satisfies_the_same_invariants() {
|
||||
let element = StandardFiniteElement::new(
|
||||
ElementType::Hex8,
|
||||
vec![
|
||||
Vector3::new(0.0, 0.0, 0.0),
|
||||
Vector3::new(2.0, 0.0, 0.0),
|
||||
Vector3::new(2.0, 1.0, 0.0),
|
||||
Vector3::new(0.0, 1.0, 0.0),
|
||||
Vector3::new(0.0, 0.0, 3.0),
|
||||
Vector3::new(2.0, 0.0, 3.0),
|
||||
Vector3::new(2.0, 1.0, 3.0),
|
||||
Vector3::new(0.0, 1.0, 3.0),
|
||||
],
|
||||
);
|
||||
let volume = 2.0 * 1.0 * 3.0;
|
||||
|
||||
let m = element.compute_element_matrices(&steel(), 0.0).unwrap();
|
||||
let k = m.stiffness_matrix;
|
||||
let mass = m.mass_matrix.expect("mass matrix required");
|
||||
|
||||
assert_eq!(k.nrows(), 8 * 3);
|
||||
|
||||
// Rigid translations.
|
||||
for d in 0..3 {
|
||||
let mut t = DVector::zeros(24);
|
||||
for node in 0..8 {
|
||||
t[node * 3 + d] = 1.0;
|
||||
}
|
||||
let f = &k * &t;
|
||||
assert!(
|
||||
f.amax() < 1e-6 * k.amax(),
|
||||
"3-D rigid translation in direction {d} produced force {:.3e}",
|
||||
f.amax()
|
||||
);
|
||||
}
|
||||
|
||||
let scale = k.amax();
|
||||
let num_zero = k
|
||||
.clone()
|
||||
.symmetric_eigenvalues()
|
||||
.iter()
|
||||
.filter(|&&e| e.abs() < 1e-9 * scale)
|
||||
.count();
|
||||
assert_eq!(
|
||||
num_zero, 6,
|
||||
"expected 6 rigid-body modes in 3-D, found {num_zero}"
|
||||
);
|
||||
|
||||
let total: f64 = mass.iter().sum();
|
||||
assert_relative_eq!(total, 3.0 * RHO * volume, max_relative = 1e-9);
|
||||
}
|
||||
@@ -58,32 +58,52 @@ mod hexahedron20_tests {
|
||||
// RED: Test Jacobian computation
|
||||
let hex20 = Hexahedron20::new();
|
||||
|
||||
// Create a regular hex with side length 2
|
||||
let mut node_coords = Vec::new();
|
||||
|
||||
// Corner nodes (8)
|
||||
for k in 0..2 {
|
||||
for j in 0..2 {
|
||||
for i in 0..2 {
|
||||
node_coords.push(Vector3::new(
|
||||
2.0 * i as f64 - 1.0,
|
||||
2.0 * j as f64 - 1.0,
|
||||
2.0 * k as f64 - 1.0,
|
||||
));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Mid-edge nodes (12) - simplified for test
|
||||
for _ in 8..20 {
|
||||
node_coords.push(Vector3::new(0.0, 0.0, 0.0));
|
||||
}
|
||||
// The reference cube spanning [-1, 1] in each direction, in the node
|
||||
// order the shape functions are written in: four bottom corners
|
||||
// counter-clockwise, four top corners, then the bottom, top and
|
||||
// vertical mid-edges.
|
||||
//
|
||||
// This previously listed the corners in lexicographic order and put
|
||||
// *all twelve* mid-edge nodes at the origin, commented "simplified
|
||||
// for test". That is not a degenerate hexahedron so much as not a
|
||||
// hexahedron: collapsing the mid-side nodes to a point makes the
|
||||
// mapping genuinely singular. It only passed because the dN/dr and
|
||||
// dN/ds derivatives carried sign errors at nodes 1, 3, 5 and 7, which
|
||||
// produced a non-zero determinant for a geometry that has none.
|
||||
let node_coords = vec![
|
||||
// Bottom corners
|
||||
Vector3::new(-1.0, -1.0, -1.0),
|
||||
Vector3::new(1.0, -1.0, -1.0),
|
||||
Vector3::new(1.0, 1.0, -1.0),
|
||||
Vector3::new(-1.0, 1.0, -1.0),
|
||||
// Top corners
|
||||
Vector3::new(-1.0, -1.0, 1.0),
|
||||
Vector3::new(1.0, -1.0, 1.0),
|
||||
Vector3::new(1.0, 1.0, 1.0),
|
||||
Vector3::new(-1.0, 1.0, 1.0),
|
||||
// Bottom mid-edges
|
||||
Vector3::new(0.0, -1.0, -1.0),
|
||||
Vector3::new(1.0, 0.0, -1.0),
|
||||
Vector3::new(0.0, 1.0, -1.0),
|
||||
Vector3::new(-1.0, 0.0, -1.0),
|
||||
// Top mid-edges
|
||||
Vector3::new(0.0, -1.0, 1.0),
|
||||
Vector3::new(1.0, 0.0, 1.0),
|
||||
Vector3::new(0.0, 1.0, 1.0),
|
||||
Vector3::new(-1.0, 0.0, 1.0),
|
||||
// Vertical mid-edges
|
||||
Vector3::new(-1.0, -1.0, 0.0),
|
||||
Vector3::new(1.0, -1.0, 0.0),
|
||||
Vector3::new(1.0, 1.0, 0.0),
|
||||
Vector3::new(-1.0, 1.0, 0.0),
|
||||
];
|
||||
|
||||
let coords = NaturalCoords::new_3d(0.0, 0.0, 0.0);
|
||||
let jac = hex20.jacobian(&coords, &node_coords).unwrap();
|
||||
|
||||
// GREEN: For regular hex, Jacobian at center should be diagonal
|
||||
assert!(jac.determinant().abs() > 1e-10); // Non-singular
|
||||
// The element occupies its own reference domain, so the mapping is the
|
||||
// identity and the Jacobian determinant is exactly 1.
|
||||
assert!((jac.determinant() - 1.0).abs() < 1e-10);
|
||||
}
|
||||
|
||||
#[test]
|
||||
|
||||
@@ -69,13 +69,29 @@ fn test_material_database() {
|
||||
|
||||
#[test]
|
||||
fn test_element_factory() {
|
||||
// Test element creation for all supported types
|
||||
// `Point`, `Line2` and `Line3` have no interpolation over an area or
|
||||
// volume, so the factory deliberately rejects them. This previously
|
||||
// required every variant of `ElementType::all()` to construct, which
|
||||
// could only pass if the factory stopped making that distinction.
|
||||
let unsupported = [ElementType::Point, ElementType::Line2, ElementType::Line3];
|
||||
|
||||
for element_type in ElementType::all() {
|
||||
let element = ElementFactory::create(element_type);
|
||||
assert!(element.is_ok(), "Failed to create {:?}", element_type);
|
||||
|
||||
let elem = element.unwrap();
|
||||
assert!(elem.num_nodes() > 0);
|
||||
if unsupported.contains(&element_type) {
|
||||
assert!(
|
||||
element.is_err(),
|
||||
"{element_type:?} has no area or volume interpolation and should be rejected"
|
||||
);
|
||||
continue;
|
||||
}
|
||||
|
||||
let elem = element.unwrap_or_else(|e| panic!("failed to create {element_type:?}: {e}"));
|
||||
assert_eq!(
|
||||
elem.num_nodes(),
|
||||
element_type.node_count(),
|
||||
"{element_type:?} reported the wrong node count"
|
||||
);
|
||||
assert!(elem.spatial_dimension() >= 1 && elem.spatial_dimension() <= 3);
|
||||
}
|
||||
}
|
||||
@@ -87,15 +103,25 @@ fn test_dof_numbering() {
|
||||
let dof_numbering =
|
||||
AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::Sequential).unwrap();
|
||||
|
||||
assert_eq!(dof_numbering.total_dofs, 27); // 9 nodes × 3 displacement components
|
||||
assert_eq!(dof_numbering.num_free_dofs(), 27);
|
||||
// 9 nodes x 2 displacement components: `generate_rectangle` builds a
|
||||
// planar mesh, which has no out-of-plane displacement to number.
|
||||
//
|
||||
// This previously expected 27, from numbering all three components on a
|
||||
// 2-D mesh. That gave every node a degree of freedom the element
|
||||
// matrices never supply, so assembly rejected every element with a
|
||||
// dimension mismatch — and `comprehensive_tdd_tests::test_dof_numbering`
|
||||
// asserts `num_nodes * 2` for the same situation, so the two tests
|
||||
// contradicted each other.
|
||||
assert_eq!(dof_numbering.total_dofs, 18);
|
||||
assert_eq!(dof_numbering.num_free_dofs(), 18);
|
||||
assert_eq!(dof_numbering.num_constrained_dofs(), 0);
|
||||
|
||||
// Test bandwidth optimized numbering
|
||||
// Test bandwidth optimized numbering. Reordering DOFs must not change how
|
||||
// many there are.
|
||||
let optimized_numbering =
|
||||
AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::BandwidthOptimized)
|
||||
.unwrap();
|
||||
assert_eq!(optimized_numbering.total_dofs, 27);
|
||||
assert_eq!(optimized_numbering.total_dofs, 18);
|
||||
}
|
||||
|
||||
#[test]
|
||||
|
||||
@@ -0,0 +1,331 @@
|
||||
//! Modal analysis against closed-form natural frequencies.
|
||||
//!
|
||||
//! This is the end-to-end check that mesh, DOF numbering, constraints,
|
||||
//! element matrices, global assembly and the eigensolver are all correct
|
||||
//! *together*. Each has its own unit tests; none of those would catch a
|
||||
//! mismatch between them, such as a mass matrix assembled in a different DOF
|
||||
//! order than its stiffness.
|
||||
//!
|
||||
//! # Why axial modes and not a cantilever
|
||||
//!
|
||||
//! The obvious benchmark is the bending frequency of a cantilever,
|
||||
//! `β₁L = 1.8751`. It is the wrong first test here. Bilinear `Quad4` elements
|
||||
//! suffer **shear locking** in bending: their assumed displacement field
|
||||
//! cannot represent pure bending without spurious shear strain, so a coarse
|
||||
//! mesh is far too stiff and reports frequencies well above the true value.
|
||||
//! A cantilever test would fail for a reason that has nothing to do with
|
||||
//! whether the code under test is correct, and tuning the tolerance until it
|
||||
//! passed would destroy its value as evidence.
|
||||
//!
|
||||
//! Longitudinal (axial) vibration has no such problem. The exact solution of
|
||||
//! the 1-D wave equation for a fixed-free bar is
|
||||
//!
|
||||
//! ```text
|
||||
//! f_n = (2n - 1) / (4L) * sqrt(E / rho), n = 1, 2, 3, ...
|
||||
//! ```
|
||||
//!
|
||||
//! and a plane-stress mesh with transverse motion suppressed reduces to
|
||||
//! exactly that problem. Linear elements with a consistent mass matrix
|
||||
//! converge to it from above at `O(h²)`, so a modest mesh lands within a
|
||||
//! fraction of a percent — tight enough that a real error cannot hide.
|
||||
//!
|
||||
//! Bending is still checked below, but as a *convergence* statement rather
|
||||
//! than a single tolerance, which is the honest way to assert on an element
|
||||
//! that is known to lock.
|
||||
|
||||
use rtx_fea::analysis::{Analysis, AnalysisConfig, AnalysisData, ModalAnalysis};
|
||||
use rtx_fea::assembly::DofComponent;
|
||||
use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, DirichletBC};
|
||||
use rtx_fea::materials::{LinearElastic, MaterialDatabase};
|
||||
use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
|
||||
|
||||
const E: f64 = 200e9;
|
||||
const RHO: f64 = 8000.0;
|
||||
const LENGTH: f64 = 1.0;
|
||||
const HEIGHT: f64 = 0.05;
|
||||
|
||||
/// A rectangular `nx` by `ny` grid of `Quad4` elements spanning
|
||||
/// `[0, LENGTH] x [0, HEIGHT]`, returned with its node grid so tests can pick
|
||||
/// out edges to constrain.
|
||||
fn bar_mesh(nx: usize, ny: usize) -> (Mesh, Vec<Vec<NodeId>>) {
|
||||
let mut mesh = Mesh::new(2).unwrap();
|
||||
|
||||
let mut grid = vec![vec![NodeId(0); ny + 1]; nx + 1];
|
||||
for (i, column) in grid.iter_mut().enumerate() {
|
||||
for (j, slot) in column.iter_mut().enumerate() {
|
||||
let x = LENGTH * i as f64 / nx as f64;
|
||||
let y = HEIGHT * j as f64 / ny as f64;
|
||||
*slot = mesh.add_node(Node::new_2d(x, y));
|
||||
}
|
||||
}
|
||||
|
||||
for i in 0..nx {
|
||||
for j in 0..ny {
|
||||
// Counter-clockwise, so the Jacobian determinant is positive.
|
||||
let nodes = vec![
|
||||
grid[i][j],
|
||||
grid[i + 1][j],
|
||||
grid[i + 1][j + 1],
|
||||
grid[i][j + 1],
|
||||
];
|
||||
let element = Element::new(ElementType::Quad4, nodes, MaterialId(0)).unwrap();
|
||||
mesh.add_element(element).unwrap();
|
||||
}
|
||||
}
|
||||
|
||||
(mesh, grid)
|
||||
}
|
||||
|
||||
/// Poisson's ratio is zero throughout.
|
||||
///
|
||||
/// This is a deliberate modelling choice, not a convenience: with `nu = 0`
|
||||
/// the axial and transverse responses decouple exactly, so the plane-stress
|
||||
/// model reduces to the 1-D bar the closed form describes. A non-zero
|
||||
/// Poisson's ratio would introduce a real physical difference between the two
|
||||
/// and the comparison would no longer be exact.
|
||||
fn steel_no_poisson() -> MaterialDatabase {
|
||||
let mut materials = MaterialDatabase::new();
|
||||
materials.add_material(
|
||||
MaterialId(0),
|
||||
LinearElastic::new(E, 0.0).with_density(RHO),
|
||||
Some("steel".to_string()),
|
||||
);
|
||||
materials
|
||||
}
|
||||
|
||||
fn frequencies_of(results: &rtx_fea::analysis::AnalysisResults) -> Vec<f64> {
|
||||
match results
|
||||
.additional_data
|
||||
.get("frequencies")
|
||||
.expect("frequencies missing")
|
||||
{
|
||||
AnalysisData::Vector(v) => v.iter().copied().collect(),
|
||||
other => panic!("frequencies had unexpected type {other:?}"),
|
||||
}
|
||||
}
|
||||
|
||||
/// Longitudinal modes of a fixed-free bar against `f_n = (2n-1)/(4L)·√(E/ρ)`.
|
||||
#[test]
|
||||
fn axial_modes_match_the_closed_form_bar() {
|
||||
let nx = 24;
|
||||
let ny = 2;
|
||||
let (mesh, grid) = bar_mesh(nx, ny);
|
||||
|
||||
let mut bcs = BoundaryConditionSet::new();
|
||||
|
||||
// Clamp the left edge axially.
|
||||
let left_edge: Vec<NodeId> = grid[0].clone();
|
||||
bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
|
||||
left_edge,
|
||||
vec![DofComponent::DisplacementX],
|
||||
0.0,
|
||||
)));
|
||||
|
||||
// Suppress transverse motion everywhere, reducing the plane-stress model
|
||||
// to the 1-D bar the closed form describes. Without this the spectrum is
|
||||
// interleaved with bending modes and the comparison is meaningless.
|
||||
let all_nodes: Vec<NodeId> = grid.iter().flatten().copied().collect();
|
||||
bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
|
||||
all_nodes,
|
||||
vec![DofComponent::DisplacementY],
|
||||
0.0,
|
||||
)));
|
||||
|
||||
let num_modes = 3;
|
||||
let mut analysis = ModalAnalysis::new(
|
||||
mesh,
|
||||
steel_no_poisson(),
|
||||
num_modes,
|
||||
AnalysisConfig::default(),
|
||||
)
|
||||
.with_boundary_conditions(bcs);
|
||||
|
||||
let results = analysis.run().expect("modal analysis failed");
|
||||
let computed = frequencies_of(&results);
|
||||
|
||||
let wave_speed = (E / RHO).sqrt();
|
||||
for n in 1..=num_modes {
|
||||
let exact = (2 * n - 1) as f64 / (4.0 * LENGTH) * wave_speed;
|
||||
let got = computed[n - 1];
|
||||
let relative_error = (got - exact).abs() / exact;
|
||||
|
||||
assert!(
|
||||
relative_error < 0.01,
|
||||
"mode {n}: computed {got:.4} Hz against exact {exact:.4} Hz \
|
||||
({:.3}% error)",
|
||||
relative_error * 100.0
|
||||
);
|
||||
|
||||
// Linear elements with a consistent mass matrix are stiffer than the
|
||||
// continuum, so the discrete frequency must come in high. Landing
|
||||
// below the exact value means something is wrong even if the
|
||||
// magnitude looks plausible.
|
||||
assert!(
|
||||
got >= exact * (1.0 - 1e-9),
|
||||
"mode {n}: computed {got:.4} Hz is below the exact {exact:.4} Hz; \
|
||||
a consistent-mass discretisation cannot be softer than the continuum"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// Refining the mesh must drive the axial error down, and at the expected
|
||||
/// second-order rate.
|
||||
///
|
||||
/// A single tolerance check can be satisfied by a wrong formula with a
|
||||
/// compensating error. A convergence *rate* cannot: it pins the
|
||||
/// discretisation itself.
|
||||
#[test]
|
||||
fn axial_frequency_converges_at_second_order() {
|
||||
let wave_speed = (E / RHO).sqrt();
|
||||
let exact = wave_speed / (4.0 * LENGTH);
|
||||
|
||||
let mut errors = Vec::new();
|
||||
for nx in [4usize, 8, 16] {
|
||||
let (mesh, grid) = bar_mesh(nx, 1);
|
||||
|
||||
let mut bcs = BoundaryConditionSet::new();
|
||||
bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
|
||||
grid[0].clone(),
|
||||
vec![DofComponent::DisplacementX],
|
||||
0.0,
|
||||
)));
|
||||
bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
|
||||
grid.iter().flatten().copied().collect(),
|
||||
vec![DofComponent::DisplacementY],
|
||||
0.0,
|
||||
)));
|
||||
|
||||
let mut analysis =
|
||||
ModalAnalysis::new(mesh, steel_no_poisson(), 1, AnalysisConfig::default())
|
||||
.with_boundary_conditions(bcs);
|
||||
let results = analysis.run().expect("modal analysis failed");
|
||||
errors.push((frequencies_of(&results)[0] - exact).abs() / exact);
|
||||
}
|
||||
|
||||
for window in errors.windows(2) {
|
||||
let rate = (window[0] / window[1]).log2();
|
||||
assert!(
|
||||
rate > 1.7,
|
||||
"halving the element size reduced the error by only 2^{rate:.2}; \
|
||||
expected close to second order. errors: {errors:?}"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// An unconstrained structure has rigid-body modes, so `K` is singular.
|
||||
///
|
||||
/// The failure must be a clear error rather than a set of near-zero
|
||||
/// eigenvalues that look like real low-frequency modes.
|
||||
#[test]
|
||||
fn unconstrained_structure_is_rejected_rather_than_silently_wrong() {
|
||||
let (mesh, _) = bar_mesh(4, 1);
|
||||
let mut analysis = ModalAnalysis::new(mesh, steel_no_poisson(), 2, AnalysisConfig::default());
|
||||
|
||||
let error = analysis
|
||||
.run()
|
||||
.expect_err("an unconstrained structure must not yield frequencies");
|
||||
let message = error.to_string().to_lowercase();
|
||||
assert!(
|
||||
message.contains("singular") || message.contains("shift"),
|
||||
"error should name the singular stiffness or the shift remedy, got: {error}"
|
||||
);
|
||||
}
|
||||
|
||||
/// Every reported natural frequency must be real and positive.
|
||||
///
|
||||
/// A constrained, positive-definite structure has no zero-frequency mode. A
|
||||
/// zero or NaN here means the constraints did not reach the assembled system
|
||||
/// or the eigenvalues came back negative.
|
||||
#[test]
|
||||
fn frequencies_are_real_and_positive() {
|
||||
let (mesh, grid) = bar_mesh(6, 2);
|
||||
|
||||
let mut bcs = BoundaryConditionSet::new();
|
||||
bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
|
||||
grid[0].clone(),
|
||||
vec![DofComponent::DisplacementX, DofComponent::DisplacementY],
|
||||
0.0,
|
||||
)));
|
||||
|
||||
let mut analysis = ModalAnalysis::new(mesh, steel_no_poisson(), 4, AnalysisConfig::default())
|
||||
.with_boundary_conditions(bcs);
|
||||
let results = analysis.run().expect("modal analysis failed");
|
||||
|
||||
let frequencies = frequencies_of(&results);
|
||||
assert_eq!(frequencies.len(), 4);
|
||||
for (i, f) in frequencies.iter().enumerate() {
|
||||
assert!(
|
||||
f.is_finite() && *f > 0.0,
|
||||
"mode {} frequency is {f}, which is not a physical frequency",
|
||||
i + 1
|
||||
);
|
||||
}
|
||||
|
||||
// Ascending, since modes are named by index.
|
||||
for pair in frequencies.windows(2) {
|
||||
assert!(
|
||||
pair[1] >= pair[0],
|
||||
"frequencies are not ascending: {frequencies:?}"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// Cantilever bending, asserted as convergence rather than as a tolerance.
|
||||
///
|
||||
/// `Quad4` locks in bending, so the coarse-mesh frequency is far too high.
|
||||
/// What must still hold is that refinement moves it monotonically *towards*
|
||||
/// the Euler-Bernoulli value `f₁ = (β₁L)²/(2πL²)·√(EI/ρA)` with
|
||||
/// `β₁L = 1.8751` — and that it approaches from above, which is the signature
|
||||
/// of locking rather than of a bug.
|
||||
#[test]
|
||||
fn cantilever_bending_converges_towards_euler_bernoulli_from_above() {
|
||||
let beta_l: f64 = 1.8751;
|
||||
// Plane stress with unit thickness: A = h, I = h³/12.
|
||||
let area = HEIGHT;
|
||||
let second_moment = HEIGHT.powi(3) / 12.0;
|
||||
let exact = beta_l.powi(2) / (2.0 * std::f64::consts::PI * LENGTH.powi(2))
|
||||
* (E * second_moment / (RHO * area)).sqrt();
|
||||
|
||||
let mut computed = Vec::new();
|
||||
for (nx, ny) in [(8usize, 2usize), (16, 4), (32, 8)] {
|
||||
let (mesh, grid) = bar_mesh(nx, ny);
|
||||
|
||||
let mut bcs = BoundaryConditionSet::new();
|
||||
bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
|
||||
grid[0].clone(),
|
||||
vec![DofComponent::DisplacementX, DofComponent::DisplacementY],
|
||||
0.0,
|
||||
)));
|
||||
|
||||
let mut analysis =
|
||||
ModalAnalysis::new(mesh, steel_no_poisson(), 1, AnalysisConfig::default())
|
||||
.with_boundary_conditions(bcs);
|
||||
let results = analysis.run().expect("modal analysis failed");
|
||||
computed.push(frequencies_of(&results)[0]);
|
||||
}
|
||||
|
||||
for (i, f) in computed.iter().enumerate() {
|
||||
assert!(
|
||||
*f > exact * 0.95,
|
||||
"mesh {i}: {f:.3} Hz is below the Euler-Bernoulli value {exact:.3} Hz; \
|
||||
a locking element cannot be softer than the beam theory it approximates"
|
||||
);
|
||||
}
|
||||
|
||||
for pair in computed.windows(2) {
|
||||
assert!(
|
||||
pair[1] <= pair[0] * 1.001,
|
||||
"refining the mesh increased the bending frequency ({:?}); \
|
||||
locking must relax under refinement, not worsen",
|
||||
computed
|
||||
);
|
||||
}
|
||||
|
||||
let coarse_error = (computed[0] - exact).abs() / exact;
|
||||
let fine_error = (computed[computed.len() - 1] - exact).abs() / exact;
|
||||
assert!(
|
||||
fine_error < coarse_error,
|
||||
"refinement did not reduce the bending error: {coarse_error:.4} -> {fine_error:.4} \
|
||||
against exact {exact:.3} Hz, computed {computed:?}"
|
||||
);
|
||||
}
|
||||
@@ -100,21 +100,23 @@ mod pyramid13_tests {
|
||||
|
||||
#[test]
|
||||
fn test_pyramid13_partition_of_unity() {
|
||||
// RED: Test that shape functions sum to 1 everywhere
|
||||
// Pyramid13 is not implemented; it reports that rather than returning
|
||||
// a basis that does not form a partition of unity.
|
||||
//
|
||||
// The previous implementation summed to 4 at the base centre, and its
|
||||
// `derivatives` allocated a 13x3 matrix then wrote rows 13 to 15 --
|
||||
// copied from a sixteen-node layout -- so this test panicked on an
|
||||
// out-of-bounds index rather than on the assertion below.
|
||||
//
|
||||
// Pyramid5 is unaffected; see `test_pyramid5_partition_of_unity`.
|
||||
let pyramid13 = Pyramid13::new();
|
||||
|
||||
let test_points = vec![
|
||||
(0.0, 0.0, 0.0), // Base center
|
||||
(0.25, 0.25, 0.25), // Quarter pyramid
|
||||
];
|
||||
|
||||
for (x, y, z) in test_points {
|
||||
for (x, y, z) in [(0.0, 0.0, 0.0), (0.25, 0.25, 0.25)] {
|
||||
let coords = NaturalCoords::new_3d(x, y, z);
|
||||
let shape = pyramid13.shape_functions(&coords).unwrap();
|
||||
|
||||
// GREEN: Sum of all shape functions should be 1
|
||||
let sum: f64 = (0..13).map(|i| shape.value(i).unwrap()).sum();
|
||||
assert!((sum - 1.0).abs() < 1e-9); // Slightly more tolerance for complex element
|
||||
let error = pyramid13
|
||||
.shape_functions(&coords)
|
||||
.expect_err("Pyramid13 must report that it is unimplemented");
|
||||
assert!(error.to_string().to_lowercase().contains("not implemented"));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -0,0 +1,157 @@
|
||||
//! Invariants every finite element shape function set must satisfy.
|
||||
//!
|
||||
//! These hold for any correct element regardless of order or geometry, so a
|
||||
//! failure localises to the basis itself rather than to an accuracy budget.
|
||||
//! Checking them across the whole element library at once is what catches an
|
||||
//! element that was written but never exercised — `Wedge15` summed to 2 at
|
||||
//! mid-height and `Quadrilateral9` emitted its functions in a node order no
|
||||
//! other element in this crate uses.
|
||||
//!
|
||||
//! Two invariants, and the second is the one that gets skipped:
|
||||
//!
|
||||
//! 1. **Partition of unity**: `Σ Nᵢ(ξ) = 1` everywhere. An interpolation that
|
||||
//! does not reproduce a constant cannot represent rigid-body motion.
|
||||
//! 2. **Vanishing derivative sum**: `Σ ∂Nᵢ/∂ξⱼ = 0` everywhere. This follows
|
||||
//! from the first by differentiation, but the derivatives are usually
|
||||
//! written out by hand separately from the values, so they can drift apart
|
||||
//! — and when they do, the shape functions are consistent while the strain
|
||||
//! they produce is not.
|
||||
|
||||
use rtx_fea::elements::shape_functions::{
|
||||
Hexahedron8, Hexahedron20, Pyramid5, Pyramid13, Quadrilateral4, Quadrilateral8, Quadrilateral9,
|
||||
ShapeFunctions, Tetrahedron4, Tetrahedron10, Triangle3, Triangle6, Wedge6, Wedge15,
|
||||
};
|
||||
|
||||
/// Sample points inside the reference domain of a quadrilateral or hexahedron,
|
||||
/// which spans `[-1, 1]` in each direction.
|
||||
const CUBE_SAMPLES: &[[f64; 3]] = &[
|
||||
[0.0, 0.0, 0.0],
|
||||
[0.5, -0.25, 0.75],
|
||||
[-0.9, 0.9, -0.3],
|
||||
[0.33, 0.67, 0.1],
|
||||
[1.0, -1.0, 1.0],
|
||||
];
|
||||
|
||||
/// Sample points inside a simplex reference domain, where the coordinates are
|
||||
/// non-negative and sum to at most one.
|
||||
const SIMPLEX_SAMPLES: &[[f64; 3]] = &[
|
||||
[0.25, 0.25, 0.25],
|
||||
[0.1, 0.2, 0.3],
|
||||
[0.0, 0.0, 0.0],
|
||||
[1.0 / 3.0, 1.0 / 3.0, 0.0],
|
||||
[0.5, 0.25, 0.125],
|
||||
];
|
||||
|
||||
/// Wedges are a triangle in `(r, s)` extruded over `t ∈ [-1, 1]`, so they need
|
||||
/// their own sample set. `t = 0` is included deliberately: it is where a
|
||||
/// missing vertical mid-edge correction shows up most strongly.
|
||||
const WEDGE_SAMPLES: &[[f64; 3]] = &[
|
||||
[1.0 / 3.0, 1.0 / 3.0, 0.0],
|
||||
[0.25, 0.25, 0.0],
|
||||
[0.2, 0.3, -0.6],
|
||||
[0.5, 0.1, 0.8],
|
||||
[0.0, 0.0, -1.0],
|
||||
];
|
||||
|
||||
fn check_element<S: ShapeFunctions>(name: &str, element: &S, samples: &[[f64; 3]]) {
|
||||
let expected_nodes = element.num_nodes();
|
||||
|
||||
for point in samples {
|
||||
let xi = &point[..];
|
||||
|
||||
let values = element
|
||||
.evaluate(xi)
|
||||
.unwrap_or_else(|e| panic!("{name}: evaluate failed at {point:?}: {e}"));
|
||||
assert_eq!(
|
||||
values.len(),
|
||||
expected_nodes,
|
||||
"{name}: returned {} values for {expected_nodes} nodes",
|
||||
values.len()
|
||||
);
|
||||
|
||||
let sum: f64 = values.iter().sum();
|
||||
assert!(
|
||||
(sum - 1.0).abs() < 1e-10,
|
||||
"{name}: shape functions sum to {sum} at {point:?}, not 1 — \
|
||||
the element cannot reproduce a constant field"
|
||||
);
|
||||
|
||||
let derivatives = element
|
||||
.derivatives(xi)
|
||||
.unwrap_or_else(|e| panic!("{name}: derivatives failed at {point:?}: {e}"));
|
||||
assert_eq!(
|
||||
derivatives.nrows(),
|
||||
expected_nodes,
|
||||
"{name}: derivative matrix has {} rows for {expected_nodes} nodes",
|
||||
derivatives.nrows()
|
||||
);
|
||||
|
||||
for direction in 0..derivatives.ncols() {
|
||||
let column_sum: f64 = derivatives.column(direction).iter().sum();
|
||||
assert!(
|
||||
column_sum.abs() < 1e-10,
|
||||
"{name}: d/d(xi_{direction}) of the shape functions sums to \
|
||||
{column_sum} at {point:?}, not 0 — the derivatives are not \
|
||||
those of the values"
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn triangles_satisfy_the_invariants() {
|
||||
check_element("Triangle3", &Triangle3, SIMPLEX_SAMPLES);
|
||||
check_element("Triangle6", &Triangle6, SIMPLEX_SAMPLES);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn quadrilaterals_satisfy_the_invariants() {
|
||||
check_element("Quadrilateral4", &Quadrilateral4, CUBE_SAMPLES);
|
||||
check_element("Quadrilateral8", &Quadrilateral8, CUBE_SAMPLES);
|
||||
check_element("Quadrilateral9", &Quadrilateral9, CUBE_SAMPLES);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn tetrahedra_satisfy_the_invariants() {
|
||||
check_element("Tetrahedron4", &Tetrahedron4, SIMPLEX_SAMPLES);
|
||||
check_element("Tetrahedron10", &Tetrahedron10, SIMPLEX_SAMPLES);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn hexahedra_satisfy_the_invariants() {
|
||||
check_element("Hexahedron8", &Hexahedron8, CUBE_SAMPLES);
|
||||
check_element("Hexahedron20", &Hexahedron20, CUBE_SAMPLES);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn wedges_satisfy_the_invariants() {
|
||||
check_element("Wedge6", &Wedge6, WEDGE_SAMPLES);
|
||||
check_element("Wedge15", &Wedge15, WEDGE_SAMPLES);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn pyramid5_satisfies_the_invariants() {
|
||||
check_element("Pyramid5", &Pyramid5, SIMPLEX_SAMPLES);
|
||||
}
|
||||
|
||||
/// `Pyramid13` is not implemented and must say so rather than return values.
|
||||
///
|
||||
/// It previously returned a basis summing to 4 at the element centre, and its
|
||||
/// derivative routine indexed past the end of the matrix it had allocated.
|
||||
/// Reporting the gap is the honest behaviour until a correct rational basis
|
||||
/// and a matching pyramid quadrature rule both exist.
|
||||
#[test]
|
||||
fn pyramid13_reports_that_it_is_unimplemented() {
|
||||
let error = Pyramid13
|
||||
.evaluate(&[0.0, 0.0, 0.0])
|
||||
.expect_err("Pyramid13 must not return shape function values");
|
||||
assert!(
|
||||
error.to_string().to_lowercase().contains("not implemented"),
|
||||
"error should say the element is unimplemented, got: {error}"
|
||||
);
|
||||
|
||||
let error = Pyramid13
|
||||
.derivatives(&[0.0, 0.0, 0.0])
|
||||
.expect_err("Pyramid13 must not return shape function derivatives");
|
||||
assert!(error.to_string().to_lowercase().contains("not implemented"));
|
||||
}
|
||||
@@ -31,21 +31,45 @@ mod standalone_tests {
|
||||
);
|
||||
}
|
||||
|
||||
/// Exercise the mesh operations rather than grepping for their names.
|
||||
///
|
||||
/// This previously searched the text of `src/mesh/mod.rs` for the strings
|
||||
/// "add_node", "add_element" and "generate_rectangle". It broke when the
|
||||
/// implementations moved into submodules, which is the smaller problem:
|
||||
/// the larger one is that a source-text search cannot distinguish a
|
||||
/// working function from one that returns zeros. Every such check in this
|
||||
/// file passed for the entire period during which element matrix
|
||||
/// computation was a stub returning `DMatrix::zeros`, quadrature returned
|
||||
/// no points at all, and cloning the material database silently dropped
|
||||
/// every material.
|
||||
#[test]
|
||||
fn test_mesh_has_real_algorithms() {
|
||||
// Verify mesh module has real implementations
|
||||
let mesh_path = include_str!("../src/mesh/mod.rs");
|
||||
use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node};
|
||||
|
||||
// Check for real mesh operations
|
||||
assert!(mesh_path.contains("add_node"), "Missing node addition");
|
||||
assert!(
|
||||
mesh_path.contains("add_element"),
|
||||
"Missing element addition"
|
||||
);
|
||||
assert!(
|
||||
mesh_path.contains("generate_rectangle"),
|
||||
"Missing mesh generation"
|
||||
);
|
||||
let mut mesh = Mesh::new(2).unwrap();
|
||||
|
||||
let n0 = mesh.add_node(Node::new_2d(0.0, 0.0));
|
||||
let n1 = mesh.add_node(Node::new_2d(1.0, 0.0));
|
||||
let n2 = mesh.add_node(Node::new_2d(1.0, 1.0));
|
||||
let n3 = mesh.add_node(Node::new_2d(0.0, 1.0));
|
||||
assert_eq!(mesh.num_nodes(), 4);
|
||||
|
||||
mesh.add_element(
|
||||
Element::new(ElementType::Quad4, vec![n0, n1, n2, n3], MaterialId(0)).unwrap(),
|
||||
)
|
||||
.unwrap();
|
||||
assert_eq!(mesh.num_elements(), 1);
|
||||
|
||||
// Nodes come back at the coordinates they went in at.
|
||||
let position = mesh.get_node(n2).expect("node 2 should exist").position();
|
||||
assert!((position.x - 1.0).abs() < 1e-12);
|
||||
assert!((position.y - 1.0).abs() < 1e-12);
|
||||
|
||||
// `generate_rectangle` takes node counts per direction, so a 3 x 3
|
||||
// grid of nodes yields 9 nodes and 2 x 2 = 4 quadrilaterals.
|
||||
let generated = Mesh::generate_rectangle(1.0, 1.0, 3, 3).unwrap();
|
||||
assert_eq!(generated.num_nodes(), 9);
|
||||
assert_eq!(generated.num_elements(), 4);
|
||||
}
|
||||
|
||||
#[test]
|
||||
|
||||
Reference in New Issue
Block a user