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