runmat_analysis_fea/assembly/
solver_solid.rs

1use std::collections::BTreeMap;
2
3use runmat_meshing_core::{ElementOrder, SolverMeshArtifact};
4
5use crate::operator::CsrMatrix;
6
7use super::{
8    elements::solid::{
9        global_stiffness_matrix as tetrahedron4_stiffness, tetrahedron10_global_stiffness_matrix,
10        SolidMaterial, Tetrahedron10ElementGeometry, Tetrahedron4ElementGeometry,
11    },
12    solid_matrix::{empty_rows, rows_to_csr, scatter_csr},
13};
14
15#[derive(Debug, Clone, Copy, PartialEq, Eq)]
16pub struct SolverSolidTopology {
17    pub dof_count: usize,
18    pub node_count: usize,
19    pub volume_element_count: usize,
20    pub order: ElementOrder,
21}
22
23#[derive(Debug, Clone, PartialEq, Eq)]
24pub enum SolverSolidAssemblyError {
25    InvalidArtifact(String),
26    UnknownElementNode { element_id: u64, node_id: u64 },
27    ElementStiffness { element_id: u64, message: String },
28}
29
30pub fn solver_solid_topology(
31    artifact: &SolverMeshArtifact,
32    base_dof_count: usize,
33) -> Result<SolverSolidTopology, SolverSolidAssemblyError> {
34    artifact
35        .validate()
36        .map_err(|failure| SolverSolidAssemblyError::InvalidArtifact(failure.to_string()))?;
37    Ok(SolverSolidTopology {
38        dof_count: artifact
39            .topology
40            .nodes
41            .len()
42            .saturating_mul(3)
43            .max(base_dof_count),
44        node_count: artifact.topology.nodes.len(),
45        volume_element_count: artifact.topology.volume_elements.len(),
46        order: artifact.resolved_request.element_order,
47    })
48}
49
50pub fn assemble_solver_solid_stiffness_csr(
51    artifact: &SolverMeshArtifact,
52    default_material: SolidMaterial,
53    materials_by_id: &BTreeMap<String, SolidMaterial>,
54    base_dof_count: usize,
55) -> Result<CsrMatrix, SolverSolidAssemblyError> {
56    let topology = solver_solid_topology(artifact, base_dof_count)?;
57    let nodes = artifact
58        .topology
59        .nodes
60        .iter()
61        .enumerate()
62        .map(|(index, node)| (node.node_id, (node.coordinates_m, index * 3)))
63        .collect::<BTreeMap<_, _>>();
64    let mut rows = empty_rows(topology.dof_count);
65    for element in &artifact.topology.volume_elements {
66        let material = materials_by_id
67            .get(&element.material_id)
68            .copied()
69            .unwrap_or(default_material);
70        match element.order {
71            ElementOrder::Tet4 => {
72                let (coordinates, offsets) =
73                    element_geometry::<4>(element.element_id, &element.node_ids, &nodes)?;
74                let stiffness = tetrahedron4_stiffness(
75                    material,
76                    Tetrahedron4ElementGeometry {
77                        nodes_m: coordinates,
78                    },
79                )
80                .map_err(|failure| element_failure(element.element_id, failure))?;
81                scatter_csr(&mut rows, &offsets, &stiffness);
82            }
83            ElementOrder::Tet10 => {
84                let (coordinates, offsets) =
85                    element_geometry::<10>(element.element_id, &element.node_ids, &nodes)?;
86                let stiffness = tetrahedron10_global_stiffness_matrix(
87                    material,
88                    Tetrahedron10ElementGeometry {
89                        nodes_m: coordinates,
90                    },
91                )
92                .map_err(|failure| element_failure(element.element_id, failure))?;
93                scatter_csr(&mut rows, &offsets, &stiffness);
94            }
95        }
96    }
97    Ok(rows_to_csr(rows))
98}
99
100type NodeIndex = BTreeMap<u64, ([f64; 3], usize)>;
101
102fn element_geometry<const N: usize>(
103    element_id: u64,
104    node_ids: &[u64],
105    nodes: &NodeIndex,
106) -> Result<([[f64; 3]; N], [usize; N]), SolverSolidAssemblyError> {
107    let mut coordinates = [[0.0; 3]; N];
108    let mut offsets = [0; N];
109    for (local, node_id) in node_ids.iter().copied().enumerate() {
110        let (point, offset) =
111            nodes
112                .get(&node_id)
113                .copied()
114                .ok_or(SolverSolidAssemblyError::UnknownElementNode {
115                    element_id,
116                    node_id,
117                })?;
118        coordinates[local] = point;
119        offsets[local] = offset;
120    }
121    Ok((coordinates, offsets))
122}
123
124fn element_failure(element_id: u64, failure: impl std::fmt::Display) -> SolverSolidAssemblyError {
125    SolverSolidAssemblyError::ElementStiffness {
126        element_id,
127        message: failure.to_string(),
128    }
129}
130
131#[cfg(test)]
132#[path = "solver_solid/tests.rs"]
133pub(crate) mod tests;
runmat_analysis_fea/assembly/solver_solid.rs

runmat_analysis_fea/assembly/
solver_solid.rs
