conspire/domain/fem/block/
mod.rs

1#[cfg(test)]
2mod test;
3
4pub mod element;
5pub mod solid;
6pub mod surface;
7pub mod thermal;
8
9use crate::{
10    fem::{
11        Elements, NodalReferenceCoordinates,
12        block::element::{
13            ElementNodalReferenceCoordinates, FiniteElement,
14            planar::PlanarElementNodalReferenceCoordinates,
15        },
16    },
17    geometry::mesh::PrimitiveConnectivity,
18    math::{
19        Quantity, Tensor, TensorRank1List, TensorRank1Vec, optimize::EqualityConstraint,
20        sparse::SparseSolver,
21    },
22    units::Volume,
23};
24use std::{
25    any::type_name,
26    fmt::{self, Debug, Formatter},
27};
28
29pub struct Block<C, F, const G: usize, const M: usize, const N: usize, const P: usize> {
30    constitutive_model: C,
31    connectivity: PrimitiveConnectivity<M, N>,
32    elements: Vec<F>,
33}
34
35impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize> Block<C, F, G, M, N, P>
36where
37    F: FiniteElement<G, M, N, P>,
38{
39    fn constitutive_model(&self) -> &C {
40        &self.constitutive_model
41    }
42    fn connectivity(&self) -> &PrimitiveConnectivity<M, N> {
43        &self.connectivity
44    }
45    fn elements(&self) -> &[F] {
46        &self.elements
47    }
48    fn element_coordinates<const D: usize, I, U>(
49        coordinates: &TensorRank1Vec<D, I, U>,
50        nodes: &[usize; N],
51    ) -> TensorRank1List<D, I, N, U> {
52        nodes
53            .iter()
54            .map(|&node| coordinates[node].clone())
55            .collect()
56    }
57    pub fn volume(&self) -> Quantity<Volume> {
58        self.elements().iter().map(|element| element.volume()).sum()
59    }
60}
61
62impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize> Debug
63    for Block<C, F, G, M, N, P>
64where
65    F: FiniteElement<G, M, N, P>,
66{
67    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
68        write!(
69            f,
70            "Block {{ constitutive model: {}, {} elements }}",
71            type_name::<C>()
72                .rsplit("::")
73                .next()
74                .unwrap()
75                .split("<")
76                .next()
77                .unwrap(),
78            self.elements().len()
79        )
80    }
81}
82
83impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize> Elements
84    for Block<C, F, G, M, N, P>
85where
86    F: FiniteElement<G, M, N, P>,
87{
88    fn node_neighbors(&self, neighbors: &mut [Vec<usize>]) {
89        add_node_neighbors(self.connectivity(), neighbors)
90    }
91}
92
93impl<C, F, const G: usize, const N: usize, const P: usize>
94    From<(
95        C,
96        PrimitiveConnectivity<3, N>,
97        &NodalReferenceCoordinates<3>,
98    )> for Block<C, F, G, 3, N, P>
99where
100    F: FiniteElement<G, 3, N, P> + From<ElementNodalReferenceCoordinates<N>>,
101{
102    fn from(
103        (constitutive_model, connectivity, coordinates): (
104            C,
105            PrimitiveConnectivity<3, N>,
106            &NodalReferenceCoordinates<3>,
107        ),
108    ) -> Self {
109        let elements = connectivity
110            .iter()
111            .map(|nodes| Self::element_coordinates(coordinates, nodes).into())
112            .collect();
113        Self {
114            constitutive_model,
115            connectivity,
116            elements,
117        }
118    }
119}
120
121impl<C, F, const G: usize, const N: usize, const P: usize>
122    From<(C, Vec<[usize; N]>, &NodalReferenceCoordinates<3>)> for Block<C, F, G, 3, N, P>
123where
124    F: FiniteElement<G, 3, N, P> + From<ElementNodalReferenceCoordinates<N>>,
125{
126    fn from(
127        (constitutive_model, connectivity, coordinates): (
128            C,
129            Vec<[usize; N]>,
130            &NodalReferenceCoordinates<3>,
131        ),
132    ) -> Self {
133        Self::from((
134            constitutive_model,
135            PrimitiveConnectivity::from(connectivity),
136            coordinates,
137        ))
138    }
139}
140
141impl<C, F, const G: usize, const N: usize, const P: usize>
142    From<(
143        C,
144        PrimitiveConnectivity<2, N>,
145        &NodalReferenceCoordinates<2>,
146    )> for Block<C, F, G, 2, N, P>
147where
148    F: FiniteElement<G, 2, N, P> + From<PlanarElementNodalReferenceCoordinates<N>>,
149{
150    fn from(
151        (constitutive_model, connectivity, coordinates): (
152            C,
153            PrimitiveConnectivity<2, N>,
154            &NodalReferenceCoordinates<2>,
155        ),
156    ) -> Self {
157        let elements = connectivity
158            .iter()
159            .map(|nodes| Self::element_coordinates(coordinates, nodes).into())
160            .collect();
161        Self {
162            constitutive_model,
163            connectivity,
164            elements,
165        }
166    }
167}
168
169impl<C, F, const G: usize, const N: usize, const P: usize>
170    From<(C, Vec<[usize; N]>, &NodalReferenceCoordinates<2>)> for Block<C, F, G, 2, N, P>
171where
172    F: FiniteElement<G, 2, N, P> + From<PlanarElementNodalReferenceCoordinates<N>>,
173{
174    fn from(
175        (constitutive_model, connectivity, coordinates): (
176            C,
177            Vec<[usize; N]>,
178            &NodalReferenceCoordinates<2>,
179        ),
180    ) -> Self {
181        Self::from((
182            constitutive_model,
183            PrimitiveConnectivity::from(connectivity),
184            coordinates,
185        ))
186    }
187}
188
189pub(crate) fn add_node_neighbors<const M: usize, const N: usize>(
190    connectivity: &PrimitiveConnectivity<M, N>,
191    neighbors: &mut [Vec<usize>],
192) {
193    connectivity.iter().for_each(|nodes| {
194        nodes.iter().for_each(|&node_a| {
195            nodes
196                .iter()
197                .for_each(|&node_b| neighbors[node_a].push(node_b))
198        })
199    })
200}
201
202pub(crate) fn finalize_node_neighbors(neighbors: &mut [Vec<usize>]) {
203    neighbors.iter_mut().for_each(|nodes| {
204        nodes.sort_unstable();
205        nodes.dedup();
206    })
207}
208
209/// The sparse solver for the positions a mesh makes nonzero.
210pub(crate) fn solver_from_neighbors(
211    neighbors: &[Vec<usize>],
212    equality_constraint: &EqualityConstraint,
213    dimension: usize,
214    symmetric: bool,
215) -> SparseSolver {
216    let number_of_nodes = neighbors.len();
217    let num_coords = dimension * number_of_nodes;
218    let mut pattern: Vec<(usize, usize)> = neighbors
219        .iter()
220        .enumerate()
221        .flat_map(|(a, nodes)| {
222            nodes.iter().flat_map(move |&b| {
223                (0..dimension).flat_map(move |i| {
224                    (0..dimension).map(move |j| (dimension * a + i, dimension * b + j))
225                })
226            })
227        })
228        .collect();
229    match equality_constraint {
230        EqualityConstraint::Fixed(indices) => {
231            let mut keep = vec![true; num_coords];
232            indices.iter().for_each(|&index| keep[index] = false);
233            let mut remap = vec![0; num_coords];
234            let mut next = 0;
235            (0..num_coords).for_each(|i| {
236                if keep[i] {
237                    remap[i] = next;
238                    next += 1;
239                }
240            });
241            pattern.retain(|&(i, j)| keep[i] && keep[j]);
242            let pattern = pattern
243                .into_iter()
244                .map(|(i, j)| (remap[i], remap[j]))
245                .collect();
246            SparseSolver::from_pattern(next, pattern, symmetric)
247        }
248        EqualityConstraint::Linear(matrix, _) => {
249            assert_eq!(matrix.width(), num_coords);
250            let num_dof = matrix.len() + matrix.width();
251            matrix.iter().enumerate().for_each(|(row, matrix_i)| {
252                let index = num_coords + row;
253                matrix_i.iter().enumerate().for_each(|(j, matrix_ij)| {
254                    if matrix_ij != &0.0 {
255                        pattern.push((index, j));
256                        pattern.push((j, index));
257                    }
258                })
259            });
260            SparseSolver::from_pattern(num_dof, pattern, symmetric)
261        }
262        EqualityConstraint::None => SparseSolver::from_pattern(num_coords, pattern, symmetric),
263    }
264}
conspire/domain/fem/block/mod.rs

conspire/domain/fem/block/
mod.rs
