runmat_analysis_fea/post/

fields.rs

use runmat_analysis_core::AnalysisField;

use crate::contracts::{
    FEA_FIELD_STRUCTURAL_BEAM_AXIAL_FORCE, FEA_FIELD_STRUCTURAL_BEAM_BENDING_MOMENT,
    FEA_FIELD_STRUCTURAL_BEAM_BENDING_STRESS, FEA_FIELD_STRUCTURAL_BEAM_SHEAR_FORCE,
    FEA_FIELD_STRUCTURAL_BEAM_TORSION_MOMENT, FEA_FIELD_STRUCTURAL_BEAM_TORSION_STRESS,
    FEA_FIELD_STRUCTURAL_DISPLACEMENT, FEA_FIELD_STRUCTURAL_EQUATION_SCALE,
    FEA_FIELD_STRUCTURAL_NODAL_VON_MISES, FEA_FIELD_STRUCTURAL_REACTION_FORCE,
    FEA_FIELD_STRUCTURAL_REACTION_MOMENT, FEA_FIELD_STRUCTURAL_RESIDUAL_NORM,
    FEA_FIELD_STRUCTURAL_ROTATION, FEA_FIELD_STRUCTURAL_SHELL_BENDING_MOMENT,
    FEA_FIELD_STRUCTURAL_SHELL_MEMBRANE_FORCE, FEA_FIELD_STRUCTURAL_SHELL_TRANSVERSE_SHEAR,
    FEA_FIELD_STRUCTURAL_SHELL_VON_MISES, FEA_FIELD_STRUCTURAL_STRAIN,
    FEA_FIELD_STRUCTURAL_STRAIN_ENERGY_DENSITY, FEA_FIELD_STRUCTURAL_STRESS,
    FEA_FIELD_STRUCTURAL_TOTAL_STRAIN_ENERGY, FEA_FIELD_STRUCTURAL_VON_MISES,
};
use crate::{
    assembly::{
        dofs::StructuralDofKind,
        elements::beam::{local_stiffness_matrix, BEAM_ELEMENT_DOF_COUNT},
        elements::solid::{strain_displacement_matrix, Tetrahedron4ElementGeometry},
        AssemblySummary, BeamRecoveryElementSummary, PrepCoordinateSummary,
        PrepRecoveryEdgeSummary, ShellRecoveryElementSummary, SolidRecoveryElementSummary,
        StructuralMaterialSummary,
    },
    operator::{apply_k, apply_k_unconstrained},
    solve::linear::LinearSolveResult,
};

const VECTOR_COMPONENT_COUNT: usize = 3;
const TENSOR_COMPONENT_COUNT: usize = 6;
type LocalTriangleCoordinates = [[f64; 2]; 3];
type LocalFrame = [[f64; 3]; 3];

pub fn recover_result_fields(
    summary: &AssemblySummary,
    solve_result: &LinearSolveResult,
) -> Vec<AnalysisField> {
    if solve_result.solution.is_empty()
        || solve_result.solution.iter().any(|value| !value.is_finite())
    {
        return empty_structural_fields();
    }

    let dof_count = summary.dof_count.max(3);
    let mut displacement_values = solve_result.solution.clone();
    if displacement_values.len() < dof_count {
        displacement_values.resize(dof_count, 0.0);
    }
    if displacement_values.is_empty() {
        displacement_values = vec![0.0; dof_count];
    }

    let node_count = structural_displacement_node_count(summary, dof_count);
    displacement_values = recover_displacement(summary, &displacement_values, node_count);

    let strain_recovery = recover_structural_strain(summary, &displacement_values);
    let element_count = strain_recovery.element_count;
    let strain_values = strain_recovery.values;
    let stress_values = recover_stress(&strain_values, summary.structural_material);
    let von_mises_values = recover_von_mises(&stress_values);
    let nodal_von_mises_values =
        recover_nodal_averaged_scalar(summary, &von_mises_values, node_count);
    let strain_energy_density_values =
        recover_strain_energy_density(&strain_values, &stress_values);
    let internal_force = apply_k_unconstrained(&summary.operator, &solve_result.solution);
    let reaction_values = recover_reaction_force(summary, &internal_force);
    let rotation_values = recover_rotation(summary, &solve_result.solution);
    let reaction_moment_values = recover_reaction_moment(summary, &internal_force);
    let strain_energy = recover_total_strain_energy(&solve_result.solution, &internal_force);
    let residual_metrics = recover_residual_metrics(summary, &solve_result.solution);

    let mut fields = vec![
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_DISPLACEMENT,
            vec![node_count, VECTOR_COMPONENT_COUNT],
            displacement_values,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_ROTATION,
            rotation_shape(summary),
            rotation_values,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_VON_MISES,
            vec![element_count],
            von_mises_values,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_STRAIN,
            vec![element_count, TENSOR_COMPONENT_COUNT],
            strain_values,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_STRESS,
            vec![element_count, TENSOR_COMPONENT_COUNT],
            stress_values,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_STRAIN_ENERGY_DENSITY,
            vec![element_count],
            strain_energy_density_values,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_NODAL_VON_MISES,
            vec![node_count],
            nodal_von_mises_values,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_REACTION_FORCE,
            reaction_shape(summary),
            reaction_values,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_REACTION_MOMENT,
            rotation_shape(summary),
            reaction_moment_values,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_TOTAL_STRAIN_ENERGY,
            vec![1],
            vec![strain_energy],
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_RESIDUAL_NORM,
            vec![1],
            vec![residual_metrics.normalized_residual_norm],
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_EQUATION_SCALE,
            vec![1],
            vec![residual_metrics.equation_scale],
        ),
    ];
    if let Some(beam_fields) = recover_beam_result_fields(summary, &solve_result.solution) {
        fields.extend(beam_fields);
    }
    if let Some(shell_fields) = recover_shell_result_fields(summary, &solve_result.solution) {
        fields.extend(shell_fields);
    }
    fields
}

pub fn recover_structural_stress_from_displacement(
    summary: &AssemblySummary,
    displacement: &[f64],
) -> Vec<f64> {
    let mut displacement_values = displacement.to_vec();
    let dof_count = summary.dof_count.max(3);
    if displacement_values.len() < dof_count {
        displacement_values.resize(dof_count, 0.0);
    }
    let node_count = dof_count.div_ceil(VECTOR_COMPONENT_COUNT).max(1);
    displacement_values.resize(node_count * VECTOR_COMPONENT_COUNT, 0.0);
    let strain_recovery = recover_structural_strain(summary, &displacement_values);
    recover_stress(&strain_recovery.values, summary.structural_material)
}

#[derive(Debug, Clone, Copy, PartialEq)]
pub struct StructuralFieldRecoveryMetrics {
    pub active_stiffness_edge_count: usize,
    pub prep_recovery_edge_count: usize,
    pub constrained_edge_count: usize,
    pub recovery_element_count: usize,
    pub solver_mesh_node_count: usize,
    pub solver_mesh_element_count: usize,
    pub max_edge_displacement_jump: f64,
    pub max_edge_strain_norm: f64,
    pub mean_edge_stiffness_ratio: f64,
    pub mean_edge_length_m: f64,
    pub strain_component_coverage_ratio: f64,
    pub element_geometry_node_count: usize,
    pub element_geometry_edge_count: usize,
    pub element_geometry_coverage_ratio: f64,
    pub basis: &'static str,
}

pub fn structural_field_recovery_metrics(
    summary: &AssemblySummary,
    displacement: &[f64],
) -> StructuralFieldRecoveryMetrics {
    let recovery_edges = structural_recovery_edges(summary);
    let strain_recovery = recover_structural_strain(summary, displacement);
    let active_stiffness_edge_count = recovery_edges.len();
    let prep_recovery_edge_count = recovery_edges
        .iter()
        .filter(|edge| edge.basis == StructuralRecoveryBasis::PrepElementConnectivity)
        .count();
    let constrained_edge_count = constrained_recovery_edge_count(summary);
    let mut max_edge_displacement_jump = 0.0_f64;
    let mut max_edge_strain_norm = 0.0_f64;
    let mut stiffness_ratio_sum = 0.0_f64;
    let mut edge_length_sum = 0.0_f64;
    for edge in &recovery_edges {
        let jump = displacement
            .get(edge.to_dof)
            .zip(displacement.get(edge.from_dof))
            .map(|(right, left)| (right - left).abs())
            .unwrap_or(0.0);
        max_edge_displacement_jump = max_edge_displacement_jump.max(jump);
        stiffness_ratio_sum += edge.stiffness_ratio;
        edge_length_sum += edge.edge_length_m;
    }
    for strain_tensor in strain_recovery.values.chunks_exact(TENSOR_COMPONENT_COUNT) {
        let strain_norm = strain_tensor
            .iter()
            .map(|value| value * value)
            .sum::<f64>()
            .sqrt();
        max_edge_strain_norm = max_edge_strain_norm.max(strain_norm);
    }
    let recovery_edge_count = recovery_edges.len();
    let mean_edge_stiffness_ratio = if recovery_edge_count == 0 {
        0.0
    } else {
        stiffness_ratio_sum / recovery_edge_count as f64
    };
    let mean_edge_length_m = if recovery_edge_count == 0 {
        0.0
    } else {
        edge_length_sum / recovery_edge_count as f64
    };
    let strain_component_coverage_ratio =
        strain_component_coverage_ratio(displacement, &recovery_edges);
    let element_geometry_node_count = summary
        .prep_coordinates
        .as_ref()
        .map(|coordinates| coordinates.element_geometry_node_count)
        .unwrap_or(0);
    let element_geometry_edge_count = summary
        .prep_coordinates
        .as_ref()
        .map(|coordinates| coordinates.element_geometry_edge_count)
        .unwrap_or(0);
    let element_geometry_coverage_ratio = summary
        .prep_coordinates
        .as_ref()
        .map(|coordinates| coordinates.element_geometry_coverage_ratio)
        .unwrap_or(0.0);
    let solver_mesh_element_count = summary.structural_solid_recovery.len();
    let solver_mesh_node_count = if solver_mesh_element_count == 0 {
        0
    } else {
        let mut nodes = summary
            .structural_solid_recovery
            .iter()
            .flat_map(|element| element.node_indices)
            .collect::<Vec<_>>();
        nodes.sort_unstable();
        nodes.dedup();
        nodes.len()
    };

    StructuralFieldRecoveryMetrics {
        active_stiffness_edge_count,
        prep_recovery_edge_count,
        constrained_edge_count,
        recovery_element_count: strain_recovery.element_count,
        solver_mesh_node_count,
        solver_mesh_element_count,
        max_edge_displacement_jump,
        max_edge_strain_norm,
        mean_edge_stiffness_ratio,
        mean_edge_length_m,
        strain_component_coverage_ratio,
        element_geometry_node_count,
        element_geometry_edge_count,
        element_geometry_coverage_ratio,
        basis: strain_recovery.basis,
    }
}

struct StructuralStrainRecovery {
    values: Vec<f64>,
    element_count: usize,
    basis: &'static str,
}

fn empty_structural_fields() -> Vec<AnalysisField> {
    vec![
        AnalysisField::host_f64(FEA_FIELD_STRUCTURAL_DISPLACEMENT, vec![0, 3], Vec::new()),
        AnalysisField::host_f64(FEA_FIELD_STRUCTURAL_ROTATION, vec![0, 3], Vec::new()),
        AnalysisField::host_f64(FEA_FIELD_STRUCTURAL_VON_MISES, vec![0], Vec::new()),
        AnalysisField::host_f64(FEA_FIELD_STRUCTURAL_STRAIN, vec![0, 6], Vec::new()),
        AnalysisField::host_f64(FEA_FIELD_STRUCTURAL_STRESS, vec![0, 6], Vec::new()),
        AnalysisField::host_f64(FEA_FIELD_STRUCTURAL_NODAL_VON_MISES, vec![0], Vec::new()),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_STRAIN_ENERGY_DENSITY,
            vec![0],
            Vec::new(),
        ),
        AnalysisField::host_f64(FEA_FIELD_STRUCTURAL_REACTION_FORCE, vec![0, 3], Vec::new()),
        AnalysisField::host_f64(FEA_FIELD_STRUCTURAL_REACTION_MOMENT, vec![0, 3], Vec::new()),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_TOTAL_STRAIN_ENERGY,
            vec![0],
            Vec::new(),
        ),
        AnalysisField::host_f64(FEA_FIELD_STRUCTURAL_RESIDUAL_NORM, vec![0], Vec::new()),
        AnalysisField::host_f64(FEA_FIELD_STRUCTURAL_EQUATION_SCALE, vec![0], Vec::new()),
    ]
}

#[derive(Debug, Clone, Copy, PartialEq)]
struct StructuralRecoveryEdge {
    from_dof: usize,
    to_dof: usize,
    component: usize,
    hop: usize,
    edge_length_m: f64,
    stiffness_ratio: f64,
    basis: StructuralRecoveryBasis,
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum StructuralRecoveryBasis {
    SolidTetrahedron4ConstantStrain,
    PrepConstantStrainBMatrix,
    PrepElementConnectivity,
    OperatorConnectivity,
}

impl StructuralRecoveryBasis {
    const fn as_str(self) -> &'static str {
        match self {
            Self::SolidTetrahedron4ConstantStrain => "solid_tetrahedron4_constant_strain",
            Self::PrepConstantStrainBMatrix => "prep_constant_strain_b_matrix",
            Self::PrepElementConnectivity => "prep_element_connectivity",
            Self::OperatorConnectivity => "operator_connectivity",
        }
    }
}

#[derive(Debug, Clone, Copy)]
struct StructuralBMatrixElement {
    nodes: [usize; 3],
    coordinates_m: [[f64; 3]; 3],
}

fn structural_recovery_edges(summary: &AssemblySummary) -> Vec<StructuralRecoveryEdge> {
    let prep_edges = prep_structural_recovery_edges(summary);
    if !prep_edges.is_empty() {
        return prep_edges;
    }

    summary
        .operator
        .stiffness_upper
        .iter()
        .enumerate()
        .filter_map(|(from_dof, stiffness)| {
            let to_dof = from_dof + 1;
            if to_dof >= summary.dof_count
                || *stiffness <= 0.0
                || summary
                    .operator
                    .constrained
                    .get(from_dof)
                    .copied()
                    .unwrap_or(false)
                || summary
                    .operator
                    .constrained
                    .get(to_dof)
                    .copied()
                    .unwrap_or(false)
            {
                return None;
            }
            let left_diag = summary
                .operator
                .stiffness_diag
                .get(from_dof)
                .copied()
                .unwrap_or(0.0);
            let right_diag = summary
                .operator
                .stiffness_diag
                .get(to_dof)
                .copied()
                .unwrap_or(0.0);
            let diag_scale = (0.5 * (left_diag + right_diag)).abs().max(1.0);
            Some(StructuralRecoveryEdge {
                from_dof,
                to_dof,
                component: from_dof % VECTOR_COMPONENT_COUNT,
                hop: 1,
                edge_length_m: 1.0,
                stiffness_ratio: (*stiffness / diag_scale).abs(),
                basis: StructuralRecoveryBasis::OperatorConnectivity,
            })
        })
        .collect()
}

fn prep_structural_recovery_edges(summary: &AssemblySummary) -> Vec<StructuralRecoveryEdge> {
    summary
        .prep_recovery_edges
        .iter()
        .filter_map(|edge| prep_recovery_edge(summary, *edge))
        .collect()
}

fn prep_recovery_edge(
    summary: &AssemblySummary,
    edge: PrepRecoveryEdgeSummary,
) -> Option<StructuralRecoveryEdge> {
    let from_dof = edge.from_dof.min(edge.to_dof);
    let to_dof = edge.from_dof.max(edge.to_dof);
    if from_dof == to_dof
        || to_dof >= summary.dof_count
        || summary
            .operator
            .constrained
            .get(from_dof)
            .copied()
            .unwrap_or(false)
        || summary
            .operator
            .constrained
            .get(to_dof)
            .copied()
            .unwrap_or(false)
    {
        return None;
    }

    let left_diag = summary
        .operator
        .stiffness_diag
        .get(from_dof)
        .copied()
        .unwrap_or(0.0);
    let right_diag = summary
        .operator
        .stiffness_diag
        .get(to_dof)
        .copied()
        .unwrap_or(0.0);
    let diag_scale = (0.5 * (left_diag + right_diag)).abs().max(1.0);
    let hop = to_dof.abs_diff(from_dof).max(1);
    let family_scale = match edge.element_family_index {
        0 => 0.95,
        1 => 1.0,
        2 => 1.05,
        3 => 1.1,
        _ => 0.9,
    };
    Some(StructuralRecoveryEdge {
        from_dof,
        to_dof,
        component: from_dof % VECTOR_COMPONENT_COUNT,
        hop,
        edge_length_m: finite_positive_or(edge.edge_length_m, hop as f64),
        stiffness_ratio: family_scale / (hop as f64 * diag_scale.sqrt().max(1.0)),
        basis: StructuralRecoveryBasis::PrepElementConnectivity,
    })
}

fn finite_positive_or(value: f64, fallback: f64) -> f64 {
    if value.is_finite() && value > 0.0 {
        value
    } else {
        fallback
    }
}

fn constrained_recovery_edge_count(summary: &AssemblySummary) -> usize {
    let prep_constrained = summary
        .prep_recovery_edges
        .iter()
        .filter(|edge| {
            let from_dof = edge.from_dof.min(edge.to_dof);
            let to_dof = edge.from_dof.max(edge.to_dof);
            to_dof < summary.dof_count
                && (summary
                    .operator
                    .constrained
                    .get(from_dof)
                    .copied()
                    .unwrap_or(false)
                    || summary
                        .operator
                        .constrained
                        .get(to_dof)
                        .copied()
                        .unwrap_or(false))
        })
        .count();
    if prep_constrained > 0 || !summary.prep_recovery_edges.is_empty() {
        return prep_constrained;
    }

    summary
        .operator
        .stiffness_upper
        .iter()
        .enumerate()
        .filter(|(from_dof, stiffness)| {
            let to_dof = from_dof + 1;
            **stiffness > 0.0
                && to_dof < summary.dof_count
                && (summary
                    .operator
                    .constrained
                    .get(*from_dof)
                    .copied()
                    .unwrap_or(false)
                    || summary
                        .operator
                        .constrained
                        .get(to_dof)
                        .copied()
                        .unwrap_or(false))
        })
        .count()
}

fn recover_structural_strain(
    summary: &AssemblySummary,
    displacement: &[f64],
) -> StructuralStrainRecovery {
    if !summary.structural_solid_recovery.is_empty() {
        return StructuralStrainRecovery {
            values: recover_solid_tetrahedron4_strain(
                displacement,
                &summary.structural_solid_recovery,
            ),
            element_count: summary.structural_solid_recovery.len().max(1),
            basis: StructuralRecoveryBasis::SolidTetrahedron4ConstantStrain.as_str(),
        };
    }

    let b_matrix_elements = prep_b_matrix_recovery_elements(summary, displacement);
    if !b_matrix_elements.is_empty() {
        return StructuralStrainRecovery {
            values: recover_b_matrix_strain(displacement, &b_matrix_elements),
            element_count: b_matrix_elements.len().max(1),
            basis: StructuralRecoveryBasis::PrepConstantStrainBMatrix.as_str(),
        };
    }

    let recovery_edges = structural_recovery_edges(summary);
    StructuralStrainRecovery {
        values: recover_edge_strain(displacement, &recovery_edges),
        element_count: recovery_edges.len().max(1),
        basis: recovery_edges
            .first()
            .map(|edge| edge.basis.as_str())
            .unwrap_or(StructuralRecoveryBasis::OperatorConnectivity.as_str()),
    }
}

fn recover_solid_tetrahedron4_strain(
    displacement: &[f64],
    elements: &[SolidRecoveryElementSummary],
) -> Vec<f64> {
    let mut strain = vec![0.0; elements.len().max(1) * TENSOR_COMPONENT_COUNT];
    for (element_index, element) in elements.iter().enumerate() {
        let Ok(b) = strain_displacement_matrix(Tetrahedron4ElementGeometry {
            nodes_m: element.coordinates_m,
        }) else {
            continue;
        };
        let mut element_displacement = [0.0_f64; 12];
        for (local_node, node_index) in element.node_indices.iter().copied().enumerate() {
            let displacement_vector = nodal_displacement(displacement, node_index);
            let base = local_node * VECTOR_COMPONENT_COUNT;
            element_displacement[base..base + VECTOR_COMPONENT_COUNT]
                .copy_from_slice(&displacement_vector);
        }
        let base = element_index * TENSOR_COMPONENT_COUNT;
        for component in 0..TENSOR_COMPONENT_COUNT {
            strain[base + component] = b[component]
                .iter()
                .zip(element_displacement.iter())
                .map(|(shape, value)| shape * value)
                .sum();
        }
    }
    strain
}

fn recover_nodal_averaged_scalar(
    summary: &AssemblySummary,
    element_values: &[f64],
    node_count: usize,
) -> Vec<f64> {
    let mut nodal = vec![0.0_f64; node_count];
    let mut counts = vec![0_usize; node_count];
    for (element_index, element) in summary.structural_solid_recovery.iter().enumerate() {
        let Some(value) = element_values.get(element_index).copied() else {
            continue;
        };
        for node_index in element.node_indices {
            if let Some(accumulator) = nodal.get_mut(node_index) {
                *accumulator += value;
                counts[node_index] += 1;
            }
        }
    }
    for (value, count) in nodal.iter_mut().zip(counts) {
        if count > 0 {
            *value /= count as f64;
        }
    }
    nodal
}

fn prep_b_matrix_recovery_elements(
    summary: &AssemblySummary,
    displacement: &[f64],
) -> Vec<StructuralBMatrixElement> {
    let Some(prep_coordinates) = summary.prep_coordinates.as_ref() else {
        return Vec::new();
    };
    if prep_coordinates.element_geometry_coverage_ratio <= 0.0
        || prep_coordinates.element_geometry_node_count < 3
        || prep_coordinates.reference_element_area_m2 <= 0.0
        || !prep_coordinates.reference_element_area_m2.is_finite()
        || !reference_coordinates_are_valid(prep_coordinates.reference_element_coordinates_m)
    {
        return Vec::new();
    }

    let node_count = displacement.len().div_ceil(VECTOR_COMPONENT_COUNT);
    let full_elements = prep_full_b_matrix_recovery_elements(summary, node_count, prep_coordinates);
    if !full_elements.is_empty() {
        return full_elements;
    }

    let sample_elements = prep_sample_b_matrix_recovery_elements(
        summary,
        node_count,
        prep_coordinates
            .element_topology_sample_element_count
            .min(4),
        prep_coordinates.element_topology_sample_edge_count.min(8),
        prep_coordinates.element_topology_sample_element_edges,
        prep_coordinates.element_topology_sample_edge_nodes,
        prep_coordinates.element_topology_sample_node_coordinates_m,
    );
    if !sample_elements.is_empty() {
        return sample_elements;
    }

    fallback_b_matrix_recovery_elements(summary, node_count, prep_coordinates)
}

fn prep_full_b_matrix_recovery_elements(
    summary: &AssemblySummary,
    node_count: usize,
    prep_coordinates: &PrepCoordinateSummary,
) -> Vec<StructuralBMatrixElement> {
    if prep_coordinates.element_topology_edge_nodes.len() < 3
        || prep_coordinates.element_topology_element_edges.is_empty()
        || prep_coordinates
            .element_topology_node_coordinates_m
            .is_empty()
    {
        return Vec::new();
    }
    prep_coordinates
        .element_topology_element_edges
        .iter()
        .filter_map(|element_edges| {
            let mut nodes = Vec::with_capacity(3);
            for edge_index in element_edges {
                let edge_index = *edge_index as usize;
                let edge_nodes = *prep_coordinates
                    .element_topology_edge_nodes
                    .get(edge_index)?;
                for node in edge_nodes {
                    let node = node as usize;
                    if node < node_count
                        && node < prep_coordinates.element_topology_node_coordinates_m.len()
                        && node_has_unconstrained_dof(summary, node)
                        && !nodes.contains(&node)
                    {
                        nodes.push(node);
                    }
                }
            }
            if nodes.len() != 3 {
                return None;
            }
            let coordinates_m = [
                prep_coordinates.element_topology_node_coordinates_m[nodes[0]],
                prep_coordinates.element_topology_node_coordinates_m[nodes[1]],
                prep_coordinates.element_topology_node_coordinates_m[nodes[2]],
            ];
            reference_coordinates_are_valid(coordinates_m).then_some(StructuralBMatrixElement {
                nodes: [nodes[0], nodes[1], nodes[2]],
                coordinates_m,
            })
        })
        .collect()
}

fn fallback_b_matrix_recovery_elements(
    summary: &AssemblySummary,
    node_count: usize,
    prep_coordinates: &PrepCoordinateSummary,
) -> Vec<StructuralBMatrixElement> {
    let mut nodes = summary
        .prep_recovery_edges
        .iter()
        .flat_map(|edge| {
            [
                edge.from_dof / VECTOR_COMPONENT_COUNT,
                edge.to_dof / VECTOR_COMPONENT_COUNT,
            ]
        })
        .filter(|node| *node < node_count && node_has_unconstrained_dof(summary, *node))
        .collect::<Vec<_>>();
    nodes.sort_unstable();
    nodes.dedup();
    if nodes.len() < 3 {
        nodes = (0..node_count)
            .filter(|node| node_has_unconstrained_dof(summary, *node))
            .take(3)
            .collect();
    }

    nodes
        .chunks_exact(3)
        .map(|chunk| StructuralBMatrixElement {
            nodes: [chunk[0], chunk[1], chunk[2]],
            coordinates_m: prep_coordinates.reference_element_coordinates_m,
        })
        .collect()
}

fn prep_sample_b_matrix_recovery_elements(
    summary: &AssemblySummary,
    node_count: usize,
    sample_element_count: usize,
    sample_edge_count: usize,
    element_edges: [[u32; 3]; 4],
    edge_nodes: [[u32; 2]; 8],
    node_coordinates_m: [[f64; 3]; 8],
) -> Vec<StructuralBMatrixElement> {
    if sample_element_count == 0 || sample_edge_count < 3 {
        return Vec::new();
    }
    element_edges
        .iter()
        .take(sample_element_count)
        .filter_map(|sample_edges| {
            let mut nodes = Vec::with_capacity(3);
            for edge_index in sample_edges {
                let edge_index = *edge_index as usize;
                if edge_index >= sample_edge_count {
                    return None;
                }
                for node in edge_nodes[edge_index] {
                    let node = node as usize;
                    if node < node_count
                        && node < node_coordinates_m.len()
                        && node_has_unconstrained_dof(summary, node)
                        && !nodes.contains(&node)
                    {
                        nodes.push(node);
                    }
                }
            }
            if nodes.len() != 3 {
                return None;
            }
            let coordinates_m = [
                node_coordinates_m[nodes[0]],
                node_coordinates_m[nodes[1]],
                node_coordinates_m[nodes[2]],
            ];
            reference_coordinates_are_valid(coordinates_m).then_some(StructuralBMatrixElement {
                nodes: [nodes[0], nodes[1], nodes[2]],
                coordinates_m,
            })
        })
        .collect()
}

fn reference_coordinates_are_valid(coordinates: [[f64; 3]; 3]) -> bool {
    coordinates.iter().flatten().all(|value| value.is_finite())
        && triangle_area_3d_m2(coordinates) > 0.0
}

fn node_has_unconstrained_dof(summary: &AssemblySummary, node: usize) -> bool {
    let base = node * VECTOR_COMPONENT_COUNT;
    (0..VECTOR_COMPONENT_COUNT).any(|component| {
        !summary
            .operator
            .constrained
            .get(base + component)
            .copied()
            .unwrap_or(false)
    })
}

fn recover_b_matrix_strain(
    displacement: &[f64],
    elements: &[StructuralBMatrixElement],
) -> Vec<f64> {
    let mut strain = vec![0.0; elements.len().max(1) * TENSOR_COMPONENT_COUNT];
    for (element_index, element) in elements.iter().enumerate() {
        let edge_strain = b_matrix_strain_tensor(displacement, element);
        let base = element_index * TENSOR_COMPONENT_COUNT;
        strain[base..base + TENSOR_COMPONENT_COUNT].copy_from_slice(&edge_strain);
    }
    strain
}

fn b_matrix_strain_tensor(displacement: &[f64], element: &StructuralBMatrixElement) -> [f64; 6] {
    let Some((local_coordinates, local_basis)) = local_triangle_coordinates(element.coordinates_m)
    else {
        return [0.0; TENSOR_COMPONENT_COUNT];
    };
    let denominator = triangle_signed_area2(local_coordinates);
    if !denominator.is_finite() || denominator.abs() <= f64::EPSILON {
        return [0.0; TENSOR_COMPONENT_COUNT];
    }

    let mut local_displacement = [[0.0_f64; 3]; 3];
    for (i, node) in element.nodes.iter().copied().enumerate() {
        let displacement_vector = nodal_displacement(displacement, node);
        local_displacement[i] = [
            dot3(displacement_vector, local_basis[0]),
            dot3(displacement_vector, local_basis[1]),
            dot3(displacement_vector, local_basis[2]),
        ];
    }

    let b = [
        local_coordinates[1][1] - local_coordinates[2][1],
        local_coordinates[2][1] - local_coordinates[0][1],
        local_coordinates[0][1] - local_coordinates[1][1],
    ];
    let c = [
        local_coordinates[2][0] - local_coordinates[1][0],
        local_coordinates[0][0] - local_coordinates[2][0],
        local_coordinates[1][0] - local_coordinates[0][0],
    ];

    let mut du_dx = 0.0_f64;
    let mut du_dy = 0.0_f64;
    let mut dv_dx = 0.0_f64;
    let mut dv_dy = 0.0_f64;
    let mut dw_dx = 0.0_f64;
    let mut dw_dy = 0.0_f64;
    for i in 0..3 {
        du_dx += b[i] * local_displacement[i][0];
        du_dy += c[i] * local_displacement[i][0];
        dv_dx += b[i] * local_displacement[i][1];
        dv_dy += c[i] * local_displacement[i][1];
        dw_dx += b[i] * local_displacement[i][2];
        dw_dy += c[i] * local_displacement[i][2];
    }

    [
        du_dx / denominator,
        dv_dy / denominator,
        0.0,
        (du_dy + dv_dx) / denominator,
        dw_dy / denominator,
        dw_dx / denominator,
    ]
}

fn local_triangle_coordinates(
    coordinates: [[f64; 3]; 3],
) -> Option<(LocalTriangleCoordinates, LocalFrame)> {
    let origin = coordinates[0];
    let edge01 = sub3(coordinates[1], origin);
    let edge02 = sub3(coordinates[2], origin);
    let e1 = normalize3(edge01)?;
    let normal = normalize3(cross3(edge01, edge02))?;
    let e2 = normalize3(cross3(normal, e1))?;
    let local = coordinates.map(|point| {
        let relative = sub3(point, origin);
        [dot3(relative, e1), dot3(relative, e2)]
    });
    Some((local, [e1, e2, normal]))
}

fn triangle_signed_area2(coordinates: [[f64; 2]; 3]) -> f64 {
    coordinates[0][0] * (coordinates[1][1] - coordinates[2][1])
        + coordinates[1][0] * (coordinates[2][1] - coordinates[0][1])
        + coordinates[2][0] * (coordinates[0][1] - coordinates[1][1])
}

fn triangle_area_3d_m2(coordinates: [[f64; 3]; 3]) -> f64 {
    let edge01 = sub3(coordinates[1], coordinates[0]);
    let edge02 = sub3(coordinates[2], coordinates[0]);
    0.5 * norm3(cross3(edge01, edge02))
}

fn sub3(left: [f64; 3], right: [f64; 3]) -> [f64; 3] {
    [left[0] - right[0], left[1] - right[1], left[2] - right[2]]
}

fn dot3(left: [f64; 3], right: [f64; 3]) -> f64 {
    left[0] * right[0] + left[1] * right[1] + left[2] * right[2]
}

fn cross3(left: [f64; 3], right: [f64; 3]) -> [f64; 3] {
    [
        left[1] * right[2] - left[2] * right[1],
        left[2] * right[0] - left[0] * right[2],
        left[0] * right[1] - left[1] * right[0],
    ]
}

fn norm3(value: [f64; 3]) -> f64 {
    dot3(value, value).sqrt()
}

fn normalize3(value: [f64; 3]) -> Option<[f64; 3]> {
    let norm = norm3(value);
    (norm.is_finite() && norm > f64::EPSILON).then_some([
        value[0] / norm,
        value[1] / norm,
        value[2] / norm,
    ])
}

fn recover_edge_strain(
    displacement: &[f64],
    recovery_edges: &[StructuralRecoveryEdge],
) -> Vec<f64> {
    let element_count = recovery_edges.len().max(1);
    let mut strain = vec![0.0; element_count * TENSOR_COMPONENT_COUNT];
    for (element_index, edge) in recovery_edges.iter().enumerate() {
        let edge_strain = edge_strain_tensor(displacement, edge);
        let base = element_index * TENSOR_COMPONENT_COUNT;
        strain[base..base + TENSOR_COMPONENT_COUNT].copy_from_slice(&edge_strain);
    }
    strain
}

fn edge_strain_tensor(displacement: &[f64], edge: &StructuralRecoveryEdge) -> [f64; 6] {
    let length = finite_positive_or(edge.edge_length_m, edge.hop.max(1) as f64);
    let from_node = edge.from_dof / VECTOR_COMPONENT_COUNT;
    let to_node = edge.to_dof / VECTOR_COMPONENT_COUNT;
    if from_node == to_node {
        let jump = displacement.get(edge.to_dof).copied().unwrap_or(0.0)
            - displacement.get(edge.from_dof).copied().unwrap_or(0.0);
        let mut tensor = [0.0; TENSOR_COMPONENT_COUNT];
        tensor[edge.component] = jump / length;
        return tensor;
    }

    let du = nodal_displacement(displacement, to_node);
    let u0 = nodal_displacement(displacement, from_node);
    let gradient = [
        (du[0] - u0[0]) / length,
        (du[1] - u0[1]) / length,
        (du[2] - u0[2]) / length,
    ];
    [
        gradient[0],
        gradient[1],
        gradient[2],
        0.5 * (gradient[0] + gradient[1]),
        0.5 * (gradient[1] + gradient[2]),
        0.5 * (gradient[0] + gradient[2]),
    ]
}

fn nodal_displacement(displacement: &[f64], node: usize) -> [f64; VECTOR_COMPONENT_COUNT] {
    let base = node * VECTOR_COMPONENT_COUNT;
    [
        displacement.get(base).copied().unwrap_or(0.0),
        displacement.get(base + 1).copied().unwrap_or(0.0),
        displacement.get(base + 2).copied().unwrap_or(0.0),
    ]
}

fn strain_component_coverage_ratio(
    displacement: &[f64],
    recovery_edges: &[StructuralRecoveryEdge],
) -> f64 {
    let expected_component_count = recovery_edges.len() * TENSOR_COMPONENT_COUNT;
    if expected_component_count == 0 {
        return 0.0;
    }
    let active_component_count = recovery_edges
        .iter()
        .map(|edge| {
            edge_strain_tensor(displacement, edge)
                .iter()
                .filter(|value| value.abs() > 0.0)
                .count()
        })
        .sum::<usize>();
    active_component_count as f64 / expected_component_count as f64
}

fn recover_stress(strain: &[f64], material: StructuralMaterialSummary) -> Vec<f64> {
    let lambda = material.lame_lambda_pa.max(0.0);
    let mu = material.shear_modulus_pa.max(0.0);
    let mut stress = vec![0.0; strain.len()];
    for (element_index, strain_tensor) in strain.chunks_exact(TENSOR_COMPONENT_COUNT).enumerate() {
        let trace = strain_tensor[0] + strain_tensor[1] + strain_tensor[2];
        let base = element_index * TENSOR_COMPONENT_COUNT;
        stress[base] = lambda * trace + 2.0 * mu * strain_tensor[0];
        stress[base + 1] = lambda * trace + 2.0 * mu * strain_tensor[1];
        stress[base + 2] = lambda * trace + 2.0 * mu * strain_tensor[2];
        stress[base + 3] = mu * strain_tensor[3];
        stress[base + 4] = mu * strain_tensor[4];
        stress[base + 5] = mu * strain_tensor[5];
    }
    stress
}

fn recover_von_mises(stress: &[f64]) -> Vec<f64> {
    stress
        .chunks_exact(TENSOR_COMPONENT_COUNT)
        .map(|tensor| {
            let sxx = tensor[0];
            let syy = tensor[1];
            let szz = tensor[2];
            let txy = tensor[3];
            let tyz = tensor[4];
            let txz = tensor[5];
            (0.5 * ((sxx - syy).powi(2) + (syy - szz).powi(2) + (szz - sxx).powi(2))
                + 3.0 * (txy.powi(2) + tyz.powi(2) + txz.powi(2)))
            .sqrt()
        })
        .collect()
}

fn recover_strain_energy_density(strain: &[f64], stress: &[f64]) -> Vec<f64> {
    strain
        .chunks_exact(TENSOR_COMPONENT_COUNT)
        .zip(stress.chunks_exact(TENSOR_COMPONENT_COUNT))
        .map(|(strain_tensor, stress_tensor)| {
            0.5 * strain_tensor
                .iter()
                .zip(stress_tensor.iter())
                .map(|(strain, stress)| strain * stress)
                .sum::<f64>()
        })
        .collect()
}

struct BeamResultValues {
    axial_force: Vec<f64>,
    shear_force: Vec<f64>,
    torsion_moment: Vec<f64>,
    bending_moment: Vec<f64>,
    bending_stress: Vec<f64>,
    torsion_stress: Vec<f64>,
}

fn recover_beam_result_fields(
    summary: &AssemblySummary,
    solution: &[f64],
) -> Option<Vec<AnalysisField>> {
    if summary.structural_beam_recovery.is_empty() {
        return None;
    }
    let values = recover_beam_result_values(summary, solution);
    let element_count = summary.structural_beam_recovery.len();
    Some(vec![
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_BEAM_AXIAL_FORCE,
            vec![element_count],
            values.axial_force,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_BEAM_SHEAR_FORCE,
            vec![element_count, 2],
            values.shear_force,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_BEAM_TORSION_MOMENT,
            vec![element_count],
            values.torsion_moment,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_BEAM_BENDING_MOMENT,
            vec![element_count, 2],
            values.bending_moment,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_BEAM_BENDING_STRESS,
            vec![element_count, 2],
            values.bending_stress,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_BEAM_TORSION_STRESS,
            vec![element_count],
            values.torsion_stress,
        ),
    ])
}

fn recover_beam_result_values(summary: &AssemblySummary, solution: &[f64]) -> BeamResultValues {
    let mut axial_force = Vec::with_capacity(summary.structural_beam_recovery.len());
    let mut shear_force = Vec::with_capacity(summary.structural_beam_recovery.len() * 2);
    let mut torsion_moment = Vec::with_capacity(summary.structural_beam_recovery.len());
    let mut bending_moment = Vec::with_capacity(summary.structural_beam_recovery.len() * 2);
    let mut bending_stress = Vec::with_capacity(summary.structural_beam_recovery.len() * 2);
    let mut torsion_stress = Vec::with_capacity(summary.structural_beam_recovery.len());

    for beam in &summary.structural_beam_recovery {
        let q_global = beam_global_displacement(summary, solution, beam);
        let q_local = multiply_beam_matrix_vector(&beam.transform_global_to_local, &q_global);
        let k_local = local_stiffness_matrix(beam.section, beam.material, beam.length_m)
            .unwrap_or([[0.0; BEAM_ELEMENT_DOF_COUNT]; BEAM_ELEMENT_DOF_COUNT]);
        let f_local = multiply_beam_matrix_vector(&k_local, &q_local);

        let axial = signed_peak_pair(f_local[0], f_local[6]);
        let shear_y = signed_peak_pair(f_local[1], f_local[7]);
        let shear_z = signed_peak_pair(f_local[2], f_local[8]);
        let torsion = signed_peak_pair(f_local[3], f_local[9]);
        let bending_y = signed_peak_pair(f_local[4], f_local[10]);
        let bending_z = signed_peak_pair(f_local[5], f_local[11]);

        axial_force.push(axial);
        shear_force.extend([shear_y, shear_z]);
        torsion_moment.push(torsion);
        bending_moment.extend([bending_y, bending_z]);
        bending_stress.extend([
            bending_stress_from_moment(bending_y, beam.section.outer_fiber_z_m, beam.section.iy_m4),
            bending_stress_from_moment(bending_z, beam.section.outer_fiber_y_m, beam.section.iz_m4),
        ]);
        torsion_stress.push(torsion_stress_from_moment(
            torsion,
            beam.section.torsion_outer_radius_m,
            beam.section.torsion_j_m4,
        ));
    }

    BeamResultValues {
        axial_force,
        shear_force,
        torsion_moment,
        bending_moment,
        bending_stress,
        torsion_stress,
    }
}

fn beam_global_displacement(
    summary: &AssemblySummary,
    solution: &[f64],
    beam: &BeamRecoveryElementSummary,
) -> [f64; BEAM_ELEMENT_DOF_COUNT] {
    let mut q = [0.0; BEAM_ELEMENT_DOF_COUNT];
    for (node_offset, node_index) in [beam.node_i_index, beam.node_j_index].iter().enumerate() {
        for (component, kind) in StructuralDofKind::ORDER.iter().enumerate() {
            let target = node_offset * StructuralDofKind::ORDER.len() + component;
            q[target] = summary
                .structural_dof_layout
                .index(*node_index, *kind)
                .and_then(|row| solution.get(row).copied())
                .unwrap_or(0.0);
        }
    }
    q
}

fn multiply_beam_matrix_vector(
    matrix: &[[f64; BEAM_ELEMENT_DOF_COUNT]; BEAM_ELEMENT_DOF_COUNT],
    vector: &[f64; BEAM_ELEMENT_DOF_COUNT],
) -> [f64; BEAM_ELEMENT_DOF_COUNT] {
    let mut out = [0.0; BEAM_ELEMENT_DOF_COUNT];
    for (row, out_value) in out.iter_mut().enumerate() {
        *out_value = matrix[row]
            .iter()
            .zip(vector.iter())
            .map(|(entry, value)| entry * value)
            .sum();
    }
    out
}

fn signed_peak_pair(a: f64, b: f64) -> f64 {
    if a.abs() >= b.abs() {
        a
    } else {
        b
    }
}

fn bending_stress_from_moment(moment_n_m: f64, outer_fiber_m: f64, inertia_m4: f64) -> f64 {
    if outer_fiber_m > 0.0 && inertia_m4 > 0.0 {
        moment_n_m * outer_fiber_m / inertia_m4
    } else {
        0.0
    }
}

fn torsion_stress_from_moment(moment_n_m: f64, outer_radius_m: f64, torsion_j_m4: f64) -> f64 {
    if outer_radius_m > 0.0 && torsion_j_m4 > 0.0 {
        moment_n_m * outer_radius_m / torsion_j_m4
    } else {
        0.0
    }
}

struct ShellResultValues {
    membrane_force: Vec<f64>,
    bending_moment: Vec<f64>,
    transverse_shear: Vec<f64>,
    von_mises: Vec<f64>,
}

fn recover_shell_result_fields(
    summary: &AssemblySummary,
    solution: &[f64],
) -> Option<Vec<AnalysisField>> {
    if summary.structural_shell_recovery.is_empty() {
        return None;
    }
    let values = recover_shell_result_values(summary, solution);
    let element_count = summary.structural_shell_recovery.len();
    Some(vec![
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_SHELL_MEMBRANE_FORCE,
            vec![element_count, 3],
            values.membrane_force,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_SHELL_BENDING_MOMENT,
            vec![element_count, 3],
            values.bending_moment,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_SHELL_TRANSVERSE_SHEAR,
            vec![element_count, 2],
            values.transverse_shear,
        ),
        AnalysisField::host_f64(
            FEA_FIELD_STRUCTURAL_SHELL_VON_MISES,
            vec![element_count],
            values.von_mises,
        ),
    ])
}

fn recover_shell_result_values(summary: &AssemblySummary, solution: &[f64]) -> ShellResultValues {
    let mut membrane_force = Vec::with_capacity(summary.structural_shell_recovery.len() * 3);
    let mut bending_moment = Vec::with_capacity(summary.structural_shell_recovery.len() * 3);
    let mut transverse_shear = Vec::with_capacity(summary.structural_shell_recovery.len() * 2);
    let mut von_mises = Vec::with_capacity(summary.structural_shell_recovery.len());

    for shell in &summary.structural_shell_recovery {
        let q = shell_node_values(summary, solution, shell);
        let edge_scale = shell_area_scale(shell.area_m2);
        let t = shell.section.thickness_m.max(1.0e-12);
        let e = shell.material.youngs_modulus_pa;
        let g = shell.material.shear_modulus_pa;
        let nu = shell.material.poisson_ratio;
        let membrane_stiffness = e * t / (1.0 - nu * nu).max(1.0e-9);
        let bending_stiffness = e * t.powi(3) / (12.0 * (1.0 - nu * nu).max(1.0e-9));
        let shear_stiffness = g * t * shell.section.shear_correction;

        let ex = component_spread(&q, 0) * edge_scale;
        let ey = component_spread(&q, 1) * edge_scale;
        let gamma_xy = (component_spread(&q, 0) + component_spread(&q, 1)) * 0.5 * edge_scale;
        let kx = component_spread(&q, 3) * edge_scale;
        let ky = component_spread(&q, 4) * edge_scale;
        let kxy = component_spread(&q, 5) * edge_scale;
        let qx =
            shear_stiffness * (component_spread(&q, 2) * edge_scale + average_component(&q, 4));
        let qy =
            shear_stiffness * (component_spread(&q, 2) * edge_scale + average_component(&q, 3));

        let nx = membrane_stiffness * (ex + nu * ey);
        let ny = membrane_stiffness * (ey + nu * ex);
        let nxy = g * t * gamma_xy;
        let mx = bending_stiffness * (kx + nu * ky);
        let my = bending_stiffness * (ky + nu * kx);
        let mxy = g * t.powi(3) * kxy / 12.0;
        let membrane_vm = (nx.powi(2) - nx * ny + ny.powi(2) + 3.0 * nxy.powi(2))
            .abs()
            .sqrt()
            / t;
        let bending_vm = 6.0
            * (mx.powi(2) - mx * my + my.powi(2) + 3.0 * mxy.powi(2))
                .abs()
                .sqrt()
            / t.powi(2);

        membrane_force.extend([nx, ny, nxy]);
        bending_moment.extend([mx, my, mxy]);
        transverse_shear.extend([qx, qy]);
        von_mises.push(membrane_vm + bending_vm);
    }

    ShellResultValues {
        membrane_force,
        bending_moment,
        transverse_shear,
        von_mises,
    }
}

fn shell_node_values(
    summary: &AssemblySummary,
    solution: &[f64],
    shell: &ShellRecoveryElementSummary,
) -> [[f64; 6]; 3] {
    let mut out = [[0.0; 6]; 3];
    for (node_offset, node_index) in shell.node_indices.iter().enumerate() {
        for (component, kind) in StructuralDofKind::ORDER.iter().enumerate() {
            out[node_offset][component] = summary
                .structural_dof_layout
                .index(*node_index, *kind)
                .and_then(|row| solution.get(row).copied())
                .unwrap_or(0.0);
        }
    }
    out
}

fn shell_area_scale(area_m2: f64) -> f64 {
    1.0 / area_m2.max(1.0e-18).sqrt()
}

fn component_spread(values: &[[f64; 6]; 3], component: usize) -> f64 {
    let min = values
        .iter()
        .map(|value| value[component])
        .fold(f64::INFINITY, f64::min);
    let max = values
        .iter()
        .map(|value| value[component])
        .fold(f64::NEG_INFINITY, f64::max);
    max - min
}

fn average_component(values: &[[f64; 6]; 3], component: usize) -> f64 {
    values.iter().map(|value| value[component]).sum::<f64>() / values.len() as f64
}

fn recover_reaction_force(summary: &AssemblySummary, internal_force: &[f64]) -> Vec<f64> {
    let mut reactions = vec![0.0; reaction_shape(summary).iter().product()];
    let mut constrained_ordinal = 0usize;
    for (dof, is_constrained) in summary.operator.constrained.iter().enumerate() {
        if !*is_constrained {
            continue;
        }
        let row = constrained_ordinal / VECTOR_COMPONENT_COUNT;
        let component = constrained_ordinal % VECTOR_COMPONENT_COUNT;
        let rhs = summary.operator.rhs.get(dof).copied().unwrap_or(0.0);
        let internal = internal_force.get(dof).copied().unwrap_or(0.0);
        reactions[row * VECTOR_COMPONENT_COUNT + component] = internal - rhs;
        constrained_ordinal += 1;
    }
    reactions
}

fn recover_rotation(summary: &AssemblySummary, solution: &[f64]) -> Vec<f64> {
    let mut rotation = vec![0.0; rotation_shape(summary).iter().product()];
    if rotation.is_empty() {
        return rotation;
    }
    for row in 0..summary.structural_dof_layout.total_dof_count() {
        let Some(address) = summary.structural_dof_layout.address(row) else {
            continue;
        };
        let Some(component) = rotational_component(address.kind) else {
            continue;
        };
        let target = address.node_index * VECTOR_COMPONENT_COUNT + component;
        if target < rotation.len() {
            rotation[target] = solution.get(row).copied().unwrap_or(0.0);
        }
    }
    rotation
}

fn recover_displacement(
    summary: &AssemblySummary,
    solution: &[f64],
    node_count: usize,
) -> Vec<f64> {
    if summary.structural_rotational_dof_count == 0 {
        let mut displacement = solution.to_vec();
        displacement.resize(node_count * VECTOR_COMPONENT_COUNT, 0.0);
        return displacement;
    }

    let mut displacement = vec![0.0; node_count * VECTOR_COMPONENT_COUNT];
    for row in 0..summary.structural_dof_layout.total_dof_count() {
        let Some(address) = summary.structural_dof_layout.address(row) else {
            continue;
        };
        let Some(component) = translational_component(address.kind) else {
            continue;
        };
        let target = address.node_index * VECTOR_COMPONENT_COUNT + component;
        if target < displacement.len() {
            displacement[target] = solution.get(row).copied().unwrap_or(0.0);
        }
    }
    displacement
}

fn translational_component(kind: StructuralDofKind) -> Option<usize> {
    match kind {
        StructuralDofKind::Ux => Some(0),
        StructuralDofKind::Uy => Some(1),
        StructuralDofKind::Uz => Some(2),
        _ => None,
    }
}

fn recover_reaction_moment(summary: &AssemblySummary, internal_force: &[f64]) -> Vec<f64> {
    let mut reactions = vec![0.0; rotation_shape(summary).iter().product()];
    if reactions.is_empty() {
        return reactions;
    }
    for (dof, is_constrained) in summary.operator.constrained.iter().enumerate() {
        if !*is_constrained {
            continue;
        }
        let Some(address) = summary.structural_dof_layout.address(dof) else {
            continue;
        };
        let Some(component) = rotational_component(address.kind) else {
            continue;
        };
        let target = address.node_index * VECTOR_COMPONENT_COUNT + component;
        if target < reactions.len() {
            let rhs = summary.operator.rhs.get(dof).copied().unwrap_or(0.0);
            let internal = internal_force.get(dof).copied().unwrap_or(0.0);
            reactions[target] = internal - rhs;
        }
    }
    reactions
}

fn rotational_component(kind: StructuralDofKind) -> Option<usize> {
    match kind {
        StructuralDofKind::Rx => Some(0),
        StructuralDofKind::Ry => Some(1),
        StructuralDofKind::Rz => Some(2),
        _ => None,
    }
}

fn rotation_shape(summary: &AssemblySummary) -> Vec<usize> {
    if summary.structural_rotational_dof_count == 0 {
        return vec![0, VECTOR_COMPONENT_COUNT];
    }
    vec![
        summary.structural_dof_layout.node_count(),
        VECTOR_COMPONENT_COUNT,
    ]
}

fn structural_displacement_node_count(summary: &AssemblySummary, dof_count: usize) -> usize {
    if summary.structural_rotational_dof_count > 0 {
        summary.structural_node_count.max(1)
    } else {
        dof_count.div_ceil(VECTOR_COMPONENT_COUNT).max(1)
    }
}

fn reaction_shape(summary: &AssemblySummary) -> Vec<usize> {
    let rows = summary
        .constrained_dof_count
        .div_ceil(VECTOR_COMPONENT_COUNT);
    vec![rows, VECTOR_COMPONENT_COUNT]
}

fn recover_total_strain_energy(displacement: &[f64], internal_force: &[f64]) -> f64 {
    0.5 * displacement
        .iter()
        .zip(internal_force.iter())
        .map(|(u, force)| u * force)
        .sum::<f64>()
}

struct StructuralResidualMetrics {
    normalized_residual_norm: f64,
    equation_scale: f64,
}

fn recover_residual_metrics(
    summary: &AssemblySummary,
    displacement: &[f64],
) -> StructuralResidualMetrics {
    let mut solution = displacement.to_vec();
    solution.resize(summary.dof_count, 0.0);
    let applied = apply_k(&summary.operator, &solution);
    let residual_norm = applied
        .iter()
        .zip(summary.operator.rhs.iter())
        .map(|(lhs, rhs)| {
            let residual = lhs - rhs;
            residual * residual
        })
        .sum::<f64>()
        .sqrt();
    let equation_scale = summary
        .operator
        .rhs
        .iter()
        .map(|value| value * value)
        .sum::<f64>()
        .sqrt()
        .max(1.0);
    StructuralResidualMetrics {
        normalized_residual_norm: residual_norm / equation_scale,
        equation_scale,
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{
        assembly::assemble_linear_system,
        fixtures::{fixture_model, FixtureId},
    };
    use runmat_meshing_core::{
        contracts::artifact::ANALYSIS_MESH_SCHEMA_VERSION, AnalysisMeshArtifact, AnalysisMeshNode,
        AnalysisMeshProvenance, AnalysisMeshQualityReport, AnalysisVolumeElement, MeshSizingField,
        VolumeElementKind,
    };

    #[test]
    fn solid_tetrahedron4_recovery_uses_solver_mesh_field_shapes() {
        let model = fixture_model(FixtureId::CantileverLinearStatic);
        let summary = assemble_linear_system(&model, None, Some(tetrahedron4_mesh()), None, None);
        let solve = LinearSolveResult {
            iterations: 1,
            residual_norm: 0.0,
            converged: true,
            host_sync_count: 0,
            solver_backend: "test".to_string(),
            device_apply_k_count: 0,
            device_apply_k_attempt_count: 0,
            solution: vec![
                0.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0, 0.02, 0.0, 0.0, 0.0, 0.03,
            ],
            solver_method: "test".to_string(),
            preconditioner: "none".to_string(),
            diagnostics: Vec::new(),
        };

        let fields = recover_result_fields(&summary, &solve);
        assert_eq!(
            field_shape(&fields, FEA_FIELD_STRUCTURAL_DISPLACEMENT),
            &[4, 3]
        );
        assert_eq!(field_shape(&fields, FEA_FIELD_STRUCTURAL_STRAIN), &[1, 6]);
        assert_eq!(field_shape(&fields, FEA_FIELD_STRUCTURAL_STRESS), &[1, 6]);
        assert_eq!(field_shape(&fields, FEA_FIELD_STRUCTURAL_VON_MISES), &[1]);
        assert_eq!(
            field_shape(&fields, FEA_FIELD_STRUCTURAL_STRAIN_ENERGY_DENSITY),
            &[1]
        );
        assert_eq!(
            field_shape(&fields, FEA_FIELD_STRUCTURAL_NODAL_VON_MISES),
            &[4]
        );
        assert!(host_field(&fields, FEA_FIELD_STRUCTURAL_STRAIN_ENERGY_DENSITY)[0] >= 0.0);
        let element_von_mises = host_field(&fields, FEA_FIELD_STRUCTURAL_VON_MISES)[0];
        assert!(host_field(&fields, FEA_FIELD_STRUCTURAL_NODAL_VON_MISES)
            .iter()
            .all(|value| (*value - element_von_mises).abs() <= 1.0e-12));
        let metrics = structural_field_recovery_metrics(&summary, &solve.solution);
        assert_eq!(metrics.basis, "solid_tetrahedron4_constant_strain");
        assert_eq!(metrics.solver_mesh_node_count, 4);
        assert_eq!(metrics.solver_mesh_element_count, 1);
    }

    #[test]
    fn finite_nonconverged_solve_still_recovers_solver_mesh_fields() {
        let model = fixture_model(FixtureId::CantileverLinearStatic);
        let summary = assemble_linear_system(&model, None, Some(tetrahedron4_mesh()), None, None);
        let solve = LinearSolveResult {
            iterations: 10,
            residual_norm: 1.0e-4,
            converged: false,
            host_sync_count: 0,
            solver_backend: "test".to_string(),
            device_apply_k_count: 0,
            device_apply_k_attempt_count: 0,
            solution: vec![
                0.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0, 0.02, 0.0, 0.0, 0.0, 0.03,
            ],
            solver_method: "test".to_string(),
            preconditioner: "none".to_string(),
            diagnostics: Vec::new(),
        };

        let fields = recover_result_fields(&summary, &solve);

        assert_eq!(field_shape(&fields, FEA_FIELD_STRUCTURAL_VON_MISES), &[1]);
        assert_eq!(field_shape(&fields, FEA_FIELD_STRUCTURAL_STRAIN), &[1, 6]);
        assert_eq!(
            field_shape(&fields, FEA_FIELD_STRUCTURAL_DISPLACEMENT),
            &[4, 3]
        );
    }

    fn field_shape<'a>(fields: &'a [AnalysisField], field_id: &str) -> &'a [usize] {
        fields
            .iter()
            .find(|field| field.field_id == field_id)
            .map(|field| field.shape.as_slice())
            .expect("field should be present")
    }

    fn host_field<'a>(fields: &'a [AnalysisField], field_id: &str) -> &'a [f64] {
        fields
            .iter()
            .find(|field| field.field_id == field_id)
            .and_then(AnalysisField::as_host_f64)
            .expect("field should be present")
    }

    fn tetrahedron4_mesh() -> AnalysisMeshArtifact {
        let mut mesh = AnalysisMeshArtifact {
            schema_version: ANALYSIS_MESH_SCHEMA_VERSION.to_string(),
            mesh_id: "unit_tetrahedron".to_string(),
            nodes: vec![
                node(1, [0.0, 0.0, 0.0]),
                node(2, [1.0, 0.0, 0.0]),
                node(3, [0.0, 1.0, 0.0]),
                node(4, [0.0, 0.0, 1.0]),
            ],
            volume_elements: vec![AnalysisVolumeElement {
                element_id: "tetrahedron_1".to_string(),
                kind: VolumeElementKind::Tetrahedron4,
                node_ids: vec![1, 2, 3, 4],
                material_region_id: "solid".to_string(),
                provenance: Vec::new(),
            }],
            boundary_faces: Vec::new(),
            boundary_edges: Vec::new(),
            quality: AnalysisMeshQualityReport::default(),
            sizing: MeshSizingField::default(),
            field_topology: Vec::new(),
            backend: Default::default(),
            adaptive_iterations: Vec::new(),
            provenance: AnalysisMeshProvenance {
                algorithm: "test".to_string(),
                source_geometry_id: "geo:test".to_string(),
                source_geometry_revision: 1,
                source_geometry_sha256: None,
            },
        };
        mesh.refresh_field_topology();
        mesh
    }

    fn node(node_id: u32, coordinates_m: [f64; 3]) -> AnalysisMeshNode {
        AnalysisMeshNode {
            node_id,
            coordinates_m,
            provenance: Vec::new(),
        }
    }
}