1#[cfg(test)]
2mod test;
3
4pub mod element;
5pub mod solid;
6pub mod surface;
7pub mod thermal;
8
9use crate::{
10 defeat_message,
11 fem::{
12 NodalReferenceCoordinates,
13 block::element::{ElementNodalReferenceCoordinates, FiniteElement},
14 },
15 math::{
16 Banded, Scalar, Tensor, TestError,
17 optimize::{
18 EqualityConstraint, FirstOrderOptimization, FirstOrderRootFinding, OptimizationError,
19 SecondOrderOptimization, ZerothOrderRootFinding,
20 },
21 },
22 mechanics::{CoordinateList, Coordinates},
23};
24use std::{
25 any::type_name,
26 fmt::{self, Debug, Display, Formatter},
27 iter::repeat_n,
28};
29
30pub type Connectivity<const N: usize> = Vec<[usize; N]>;
31
32pub struct Block<C, F, const G: usize, const M: usize, const N: usize, const P: usize> {
33 constitutive_model: C,
34 connectivity: Connectivity<N>,
35 coordinates: NodalReferenceCoordinates,
36 elements: Vec<F>,
37}
38
39impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize> Block<C, F, G, M, N, P>
40where
41 F: FiniteElement<G, M, N, P>,
42{
43 fn constitutive_model(&self) -> &C {
44 &self.constitutive_model
45 }
46 fn connectivity(&self) -> &Connectivity<N> {
47 &self.connectivity
48 }
49 fn coordinates(&self) -> &NodalReferenceCoordinates {
50 &self.coordinates
51 }
52 fn elements(&self) -> &[F] {
53 &self.elements
54 }
55 fn element_coordinates<const I: usize>(
56 coordinates: &Coordinates<I>,
57 nodes: &[usize; N],
58 ) -> CoordinateList<I, N> {
59 nodes
60 .iter()
61 .map(|&node| coordinates[node].clone())
62 .collect()
63 }
64 pub fn volume(&self) -> Scalar {
65 self.elements().iter().map(|element| element.volume()).sum()
66 }
67}
68
69impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize> Debug
70 for Block<C, F, G, M, N, P>
71where
72 F: FiniteElement<G, M, N, P>,
73{
74 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
75 write!(
76 f,
77 "Block {{ constitutive model: {}, {} elements }}",
78 type_name::<C>()
79 .rsplit("::")
80 .next()
81 .unwrap()
82 .split("<")
83 .next()
84 .unwrap(),
85 self.elements().len()
86 )
87 }
88}
89
90pub trait FiniteElementBlock<C, F, const G: usize, const N: usize>
91where
92 Self: From<(C, Connectivity<N>, NodalReferenceCoordinates)>,
93{
94 fn reset(&mut self);
95}
96
97impl<C, F, const G: usize, const N: usize, const P: usize>
98 From<(C, Connectivity<N>, NodalReferenceCoordinates)> 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 Connectivity<N>,
106 NodalReferenceCoordinates,
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 coordinates,
117 elements,
118 }
119 }
120}
121
122impl<C, F, const G: usize, const N: usize, const P: usize> FiniteElementBlock<C, F, G, N>
123 for Block<C, F, G, 3, N, P>
124where
125 F: Default + FiniteElement<G, 3, N, P> + From<ElementNodalReferenceCoordinates<N>>,
126{
127 fn reset(&mut self) {
128 self.elements
129 .iter_mut()
130 .for_each(|element| *element = F::default())
131 }
132}
133
134pub enum FiniteElementBlockError {
135 Upstream(String, String),
136}
137
138impl From<FiniteElementBlockError> for String {
139 fn from(error: FiniteElementBlockError) -> Self {
140 match error {
141 FiniteElementBlockError::Upstream(error, block) => {
142 format!(
143 "{error}\x1b[0;91m\n\
144 In finite element block: {block}."
145 )
146 }
147 }
148 }
149}
150
151impl From<FiniteElementBlockError> for TestError {
152 fn from(error: FiniteElementBlockError) -> Self {
153 Self {
154 message: error.to_string(),
155 }
156 }
157}
158
159impl Debug for FiniteElementBlockError {
160 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
161 let error = match self {
162 Self::Upstream(error, block) => {
163 format!(
164 "{error}\x1b[0;91m\n\
165 In block: {block}."
166 )
167 }
168 };
169 write!(f, "\n{error}\n\x1b[0;2;31m{}\x1b[0m\n", defeat_message())
170 }
171}
172
173impl Display for FiniteElementBlockError {
174 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
175 let error = match self {
176 Self::Upstream(error, block) => {
177 format!(
178 "{error}\x1b[0;91m\n\
179 In block: {block}."
180 )
181 }
182 };
183 write!(f, "{error}\x1b[0m")
184 }
185}
186
187pub trait ZerothOrderRoot<C, E, const G: usize, const M: usize, const N: usize, X> {
188 fn root(
189 &self,
190 equality_constraint: EqualityConstraint,
191 solver: impl ZerothOrderRootFinding<X>,
192 ) -> Result<X, OptimizationError>;
193}
194
195pub trait FirstOrderRoot<C, E, const G: usize, const M: usize, const N: usize, F, J, X> {
196 fn root(
197 &self,
198 equality_constraint: EqualityConstraint,
199 solver: impl FirstOrderRootFinding<F, J, X>,
200 ) -> Result<X, OptimizationError>;
201}
202
203pub trait FirstOrderMinimize<C, E, const G: usize, const M: usize, const N: usize, X> {
204 fn minimize(
205 &self,
206 equality_constraint: EqualityConstraint,
207 solver: impl FirstOrderOptimization<Scalar, X>,
208 ) -> Result<X, OptimizationError>;
209}
210
211pub trait SecondOrderMinimize<C, E, const G: usize, const M: usize, const N: usize, J, H, X> {
212 fn minimize(
213 &self,
214 equality_constraint: EqualityConstraint,
215 solver: impl SecondOrderOptimization<Scalar, J, H, X>,
216 ) -> Result<X, OptimizationError>;
217}
218
219fn band<const N: usize>(
220 connectivity: &Connectivity<N>,
221 equality_constraint: &EqualityConstraint,
222 number_of_nodes: usize,
223 dimension: usize,
224) -> Banded {
225 match equality_constraint {
226 EqualityConstraint::Fixed(indices) => {
227 let neighbors: Vec<Vec<usize>> = invert(connectivity, number_of_nodes)
228 .iter()
229 .map(|elements| {
230 let mut nodes: Vec<usize> = elements
231 .iter()
232 .flat_map(|&element| connectivity[element])
233 .collect();
234 nodes.sort();
235 nodes.dedup();
236 nodes
237 })
238 .collect();
239 let structure: Vec<Vec<bool>> = neighbors
240 .iter()
241 .map(|nodes| (0..number_of_nodes).map(|b| nodes.contains(&b)).collect())
242 .collect();
243 let structure_nd: Vec<Vec<bool>> = structure
244 .iter()
245 .flat_map(|row| {
246 repeat_n(
247 row.iter()
248 .flat_map(|entry| repeat_n(*entry, dimension))
249 .collect(),
250 dimension,
251 )
252 })
253 .collect();
254 let mut keep = vec![true; structure_nd.len()];
255 indices.iter().for_each(|&index| keep[index] = false);
256 let banded = structure_nd
257 .into_iter()
258 .zip(keep.iter())
259 .filter(|(_, keep)| **keep)
260 .map(|(structure_nd_a, _)| {
261 structure_nd_a
262 .into_iter()
263 .zip(keep.iter())
264 .filter(|(_, keep)| **keep)
265 .map(|(structure_nd_ab, _)| structure_nd_ab)
266 .collect::<Vec<bool>>()
267 })
268 .collect::<Vec<Vec<bool>>>();
269 Banded::from(banded)
270 }
271 EqualityConstraint::Linear(matrix, _) => {
272 let neighbors: Vec<Vec<usize>> = invert(connectivity, number_of_nodes)
273 .iter()
274 .map(|elements| {
275 let mut nodes: Vec<usize> = elements
276 .iter()
277 .flat_map(|&element| connectivity[element])
278 .collect();
279 nodes.sort();
280 nodes.dedup();
281 nodes
282 })
283 .collect();
284 let structure: Vec<Vec<bool>> = neighbors
285 .iter()
286 .map(|nodes| (0..number_of_nodes).map(|b| nodes.contains(&b)).collect())
287 .collect();
288 let structure_nd: Vec<Vec<bool>> = structure
289 .iter()
290 .flat_map(|row| {
291 repeat_n(
292 row.iter()
293 .flat_map(|entry| repeat_n(*entry, dimension))
294 .collect(),
295 dimension,
296 )
297 })
298 .collect();
299 let num_coords = dimension * number_of_nodes;
300 assert_eq!(matrix.width(), num_coords);
301 let num_dof = matrix.len() + matrix.width();
302 let mut banded = vec![vec![false; num_dof]; num_dof];
303 structure_nd
304 .iter()
305 .zip(banded.iter_mut())
306 .for_each(|(structure_nd_i, banded_i)| {
307 structure_nd_i
308 .iter()
309 .zip(banded_i.iter_mut())
310 .for_each(|(structure_nd_ij, banded_ij)| *banded_ij = *structure_nd_ij)
311 });
312 let mut index = num_coords;
313 matrix.iter().for_each(|matrix_i| {
314 matrix_i.iter().enumerate().for_each(|(j, matrix_ij)| {
315 if matrix_ij != &0.0 {
316 banded[index][j] = true;
317 banded[j][index] = true;
318 index += 1;
319 }
320 })
321 });
322 Banded::from(banded)
323 }
324 EqualityConstraint::None => {
325 let neighbors: Vec<Vec<usize>> = invert(connectivity, number_of_nodes)
326 .iter()
327 .map(|elements| {
328 let mut nodes: Vec<usize> = elements
329 .iter()
330 .flat_map(|&element| connectivity[element])
331 .collect();
332 nodes.sort();
333 nodes.dedup();
334 nodes
335 })
336 .collect();
337 let structure: Vec<Vec<bool>> = neighbors
338 .iter()
339 .map(|nodes| (0..number_of_nodes).map(|b| nodes.contains(&b)).collect())
340 .collect();
341 let structure_nd: Vec<Vec<bool>> = structure
342 .iter()
343 .flat_map(|row| {
344 repeat_n(
345 row.iter()
346 .flat_map(|entry| repeat_n(*entry, dimension))
347 .collect(),
348 dimension,
349 )
350 })
351 .collect();
352 Banded::from(structure_nd)
353 }
354 }
355}
356
357fn invert<const N: usize>(
358 connectivity: &Connectivity<N>,
359 number_of_nodes: usize,
360) -> Vec<Vec<usize>> {
361 let mut inverse_connectivity = vec![vec![]; number_of_nodes];
362 connectivity
363 .iter()
364 .enumerate()
365 .for_each(|(element, nodes)| {
366 nodes
367 .iter()
368 .for_each(|&node| inverse_connectivity[node].push(element))
369 });
370 inverse_connectivity
371}