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