conspire/domain/partition/feti/block/element/
mod.rs

1use crate::{
2    domain::{Blocks, ElementModelError, Model, NodalCoordinates},
3    geometry::mesh::Partition,
4    math::{Scalar, SquareMatrix, Tensor, Vector},
5};
6use std::{array::from_fn, collections::HashMap};
7
8/// The tangent of a Newton step left unassembled.
9///
10/// Each element's stiffness and force, in the order the elements were meshed.
11/// A decomposition assigns whole elements to subdomains, so this is what one
12/// needs to build each subdomain's own system.
13pub struct ElementSystems {
14    pub(crate) positions: Vec<[f64; 3]>,
15    pub(crate) elements: Vec<ElementSystem>,
16}
17
18pub(crate) struct ElementSystem {
19    pub(crate) nodes: Vec<usize>,
20    pub(crate) stiffness: SquareMatrix,
21    pub(crate) force: Vector,
22}
23
24impl ElementSystem {
25    pub(crate) fn pack(
26        nodes: Vec<usize>,
27        force: impl Fn(usize, usize) -> Scalar,
28        stiffness: impl Fn(usize, usize, usize, usize) -> Scalar,
29    ) -> Self {
30        const D: usize = 3;
31        let number_of_nodes = nodes.len();
32        let mut packed_stiffness = SquareMatrix::zero(D * number_of_nodes);
33        let mut packed_force = Vector::zero(D * number_of_nodes);
34        (0..number_of_nodes).for_each(|a| {
35            (0..D).for_each(|i| packed_force[D * a + i] = force(a, i));
36            (0..number_of_nodes).for_each(|b| {
37                (0..D).for_each(|i| {
38                    (0..D).for_each(|j| {
39                        packed_stiffness[D * a + i][D * b + j] = stiffness(a, b, i, j)
40                    })
41                })
42            })
43        });
44        Self {
45            nodes,
46            stiffness: packed_stiffness,
47            force: packed_force,
48        }
49    }
50}
51
52pub(crate) fn positions<const D: usize>(nodal_coordinates: &NodalCoordinates<D>) -> Vec<[f64; D]> {
53    nodal_coordinates
54        .iter()
55        .map(|coordinate| from_fn(|axis| coordinate[axis].value()))
56        .collect()
57}
58
59/// Elements that can hand out their systems one by one.
60///
61/// This is what a decomposed solve asks of a model, and a model that cannot
62/// answer it, having couplings that cross any cut, is refused at compile time.
63pub trait DecomposableElements {
64    fn element_systems(
65        &self,
66        nodal_coordinates: &NodalCoordinates<3>,
67    ) -> Result<ElementSystems, ElementModelError>;
68}
69
70impl<B> DecomposableElements for Model<B, 3>
71where
72    B: DecomposableElements,
73{
74    fn element_systems(
75        &self,
76        nodal_coordinates: &NodalCoordinates<3>,
77    ) -> Result<ElementSystems, ElementModelError> {
78        self.blocks.element_systems(nodal_coordinates)
79    }
80}
81
82impl<B1, B2> DecomposableElements for Blocks<B1, B2>
83where
84    B1: DecomposableElements,
85    B2: DecomposableElements,
86{
87    fn element_systems(
88        &self,
89        nodal_coordinates: &NodalCoordinates<3>,
90    ) -> Result<ElementSystems, ElementModelError> {
91        let mut systems = self.0.element_systems(nodal_coordinates)?;
92        systems
93            .elements
94            .extend(self.1.element_systems(nodal_coordinates)?.elements);
95        Ok(systems)
96    }
97}
98
99impl ElementSystems {
100    pub(crate) fn positions(&self) -> &[[f64; 3]] {
101        &self.positions
102    }
103    #[allow(clippy::type_complexity)]
104    pub(crate) fn subdomains(
105        &self,
106        partition: &Partition,
107    ) -> Result<(Vec<SquareMatrix>, Vec<Vector>), String> {
108        const D: usize = 3;
109        if partition.elements_parts().len() != self.elements.len() {
110            return Err("The partition must assign every element to a subdomain.".to_string());
111        }
112        Ok((0..partition.number_of_parts())
113            .map(|part| {
114                let nodes = partition.part_nodes(part);
115                let local: HashMap<usize, usize> = nodes
116                    .iter()
117                    .enumerate()
118                    .map(|(local, &node)| (node, local))
119                    .collect();
120                let mut stiffness = SquareMatrix::zero(D * nodes.len());
121                let mut force = Vector::zero(D * nodes.len());
122                partition.part_elements(part).iter().for_each(|&element| {
123                    let element = &self.elements[element];
124                    let dofs: Vec<usize> = element
125                        .nodes
126                        .iter()
127                        .flat_map(|node| {
128                            let base = D * local[node];
129                            (0..D).map(move |i| base + i)
130                        })
131                        .collect();
132                    dofs.iter().enumerate().for_each(|(a, &dof_a)| {
133                        force[dof_a] += element.force[a];
134                        dofs.iter().enumerate().for_each(|(b, &dof_b)| {
135                            stiffness[dof_a][dof_b] += element.stiffness[a][b]
136                        })
137                    })
138                });
139                (stiffness, force)
140            })
141            .unzip())
142    }
143}
conspire/domain/partition/feti/block/element/mod.rs

conspire/domain/partition/feti/block/element/
mod.rs
