conspire/domain/fem/block/thermal/conduction/
mod.rs

1#[cfg(test)]
2pub mod test;
3
4use crate::{
5    constitutive::thermal::conduction::ThermalConduction,
6    fem::block::{
7        Block, FiniteElementBlockError, FirstOrderMinimize, FirstOrderRoot, SecondOrderMinimize,
8        ZerothOrderRoot, band,
9        element::{FiniteElementError, thermal::conduction::ThermalConductionFiniteElement},
10        thermal::{NodalTemperatures, ThermalFiniteElementBlock},
11    },
12    math::{
13        Scalar, SquareMatrix, Tensor, Vector,
14        optimize::{
15            EqualityConstraint, FirstOrderOptimization, FirstOrderRootFinding, OptimizationError,
16            SecondOrderOptimization, ZerothOrderRootFinding,
17        },
18    },
19};
20
21pub type NodalForcesThermal = Vector;
22pub type NodalStiffnessesThermal = SquareMatrix;
23
24pub trait ThermalConductionFiniteElementBlock<
25    C,
26    F,
27    const G: usize,
28    const M: usize,
29    const N: usize,
30    const P: usize,
31> where
32    C: ThermalConduction,
33    F: ThermalConductionFiniteElement<C, G, M, N, P>,
34{
35    fn potential(
36        &self,
37        nodal_temperatures: &NodalTemperatures,
38    ) -> Result<Scalar, FiniteElementBlockError>;
39    fn nodal_forces(
40        &self,
41        nodal_temperatures: &NodalTemperatures,
42    ) -> Result<NodalForcesThermal, FiniteElementBlockError>;
43    fn nodal_stiffnesses(
44        &self,
45        nodal_temperatures: &NodalTemperatures,
46    ) -> Result<NodalStiffnessesThermal, FiniteElementBlockError>;
47}
48
49impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
50    ThermalConductionFiniteElementBlock<C, F, G, M, N, P> for Block<C, F, G, M, N, P>
51where
52    C: ThermalConduction,
53    F: ThermalConductionFiniteElement<C, G, M, N, P>,
54{
55    fn potential(
56        &self,
57        nodal_temperatures: &NodalTemperatures,
58    ) -> Result<Scalar, FiniteElementBlockError> {
59        match self
60            .elements()
61            .iter()
62            .zip(self.connectivity())
63            .map(|(element, element_connectivity)| {
64                element.potential(
65                    self.constitutive_model(),
66                    &self.nodal_temperatures_element(element_connectivity, nodal_temperatures),
67                )
68            })
69            .sum()
70        {
71            Ok(potential) => Ok(potential),
72            Err(error) => Err(FiniteElementBlockError::Upstream(
73                format!("{error}"),
74                format!("{self:?}"),
75            )),
76        }
77    }
78    fn nodal_forces(
79        &self,
80        nodal_temperatures: &NodalTemperatures,
81    ) -> Result<NodalForcesThermal, FiniteElementBlockError> {
82        let mut nodal_forces = NodalForcesThermal::zero(nodal_temperatures.len());
83        match self
84            .elements()
85            .iter()
86            .zip(self.connectivity())
87            .try_for_each(|(element, element_connectivity)| {
88                element
89                    .nodal_forces(
90                        self.constitutive_model(),
91                        &self.nodal_temperatures_element(element_connectivity, nodal_temperatures),
92                    )?
93                    .iter()
94                    .zip(element_connectivity)
95                    .for_each(|(nodal_force, &node)| nodal_forces[node] += nodal_force);
96                Ok::<(), FiniteElementError>(())
97            }) {
98            Ok(()) => Ok(nodal_forces),
99            Err(error) => Err(FiniteElementBlockError::Upstream(
100                format!("{error}"),
101                format!("{self:?}"),
102            )),
103        }
104    }
105    fn nodal_stiffnesses(
106        &self,
107        nodal_temperatures: &NodalTemperatures,
108    ) -> Result<NodalStiffnessesThermal, FiniteElementBlockError> {
109        let mut nodal_stiffnesses = NodalStiffnessesThermal::zero(nodal_temperatures.len());
110        match self
111            .elements()
112            .iter()
113            .zip(self.connectivity())
114            .try_for_each(|(element, element_connectivity)| {
115                element
116                    .nodal_stiffnesses(
117                        self.constitutive_model(),
118                        &self.nodal_temperatures_element(element_connectivity, nodal_temperatures),
119                    )?
120                    .iter()
121                    .zip(element_connectivity)
122                    .for_each(|(object, &node_a)| {
123                        object.iter().zip(element_connectivity).for_each(
124                            |(nodal_stiffness, &node_b)| {
125                                nodal_stiffnesses[node_a][node_b] += nodal_stiffness
126                            },
127                        )
128                    });
129                Ok::<(), FiniteElementError>(())
130            }) {
131            Ok(()) => Ok(nodal_stiffnesses),
132            Err(error) => Err(FiniteElementBlockError::Upstream(
133                format!("{error}"),
134                format!("{self:?}"),
135            )),
136        }
137    }
138}
139
140impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
141    ZerothOrderRoot<C, F, G, M, N, NodalTemperatures> for Block<C, F, G, M, N, P>
142where
143    C: ThermalConduction,
144    F: ThermalConductionFiniteElement<C, G, M, N, P>,
145{
146    fn root(
147        &self,
148        equality_constraint: EqualityConstraint,
149        solver: impl ZerothOrderRootFinding<NodalTemperatures>,
150    ) -> Result<NodalTemperatures, OptimizationError> {
151        solver.root(
152            |nodal_temperatures: &NodalTemperatures| Ok(self.nodal_forces(nodal_temperatures)?),
153            NodalTemperatures::zero(self.coordinates().len()),
154            equality_constraint,
155        )
156    }
157}
158
159impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
160    FirstOrderRoot<C, F, G, M, N, NodalForcesThermal, NodalStiffnessesThermal, NodalTemperatures>
161    for Block<C, F, G, M, N, P>
162where
163    C: ThermalConduction,
164    F: ThermalConductionFiniteElement<C, G, M, N, P>,
165{
166    fn root(
167        &self,
168        equality_constraint: EqualityConstraint,
169        solver: impl FirstOrderRootFinding<
170            NodalForcesThermal,
171            NodalStiffnessesThermal,
172            NodalTemperatures,
173        >,
174    ) -> Result<NodalTemperatures, OptimizationError> {
175        solver.root(
176            |nodal_temperatures: &NodalTemperatures| Ok(self.nodal_forces(nodal_temperatures)?),
177            |nodal_temperatures: &NodalTemperatures| {
178                Ok(self.nodal_stiffnesses(nodal_temperatures)?)
179            },
180            NodalTemperatures::zero(self.coordinates().len()),
181            equality_constraint,
182        )
183    }
184}
185
186impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
187    FirstOrderMinimize<C, F, G, M, N, NodalTemperatures> for Block<C, F, G, M, N, P>
188where
189    C: ThermalConduction,
190    F: ThermalConductionFiniteElement<C, G, M, N, P>,
191{
192    fn minimize(
193        &self,
194        equality_constraint: EqualityConstraint,
195        solver: impl FirstOrderOptimization<Scalar, NodalTemperatures>,
196    ) -> Result<NodalTemperatures, OptimizationError> {
197        solver.minimize(
198            |nodal_temperatures: &NodalTemperatures| Ok(self.potential(nodal_temperatures)?),
199            |nodal_temperatures: &NodalTemperatures| Ok(self.nodal_forces(nodal_temperatures)?),
200            NodalTemperatures::zero(self.coordinates().len()),
201            equality_constraint,
202        )
203    }
204}
205
206impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
207    SecondOrderMinimize<
208        C,
209        F,
210        G,
211        M,
212        N,
213        NodalForcesThermal,
214        NodalStiffnessesThermal,
215        NodalTemperatures,
216    > for Block<C, F, G, M, N, P>
217where
218    C: ThermalConduction,
219    F: ThermalConductionFiniteElement<C, G, M, N, P>,
220{
221    fn minimize(
222        &self,
223        equality_constraint: EqualityConstraint,
224        solver: impl SecondOrderOptimization<
225            Scalar,
226            NodalForcesThermal,
227            NodalStiffnessesThermal,
228            NodalTemperatures,
229        >,
230    ) -> Result<NodalTemperatures, OptimizationError> {
231        let banded = band(
232            self.connectivity(),
233            &equality_constraint,
234            self.coordinates().len(),
235            1,
236        );
237        solver.minimize(
238            |nodal_temperatures: &NodalTemperatures| Ok(self.potential(nodal_temperatures)?),
239            |nodal_temperatures: &NodalTemperatures| Ok(self.nodal_forces(nodal_temperatures)?),
240            |nodal_temperatures: &NodalTemperatures| {
241                Ok(self.nodal_stiffnesses(nodal_temperatures)?)
242            },
243            NodalTemperatures::zero(self.coordinates().len()),
244            equality_constraint,
245            Some(banded),
246        )
247    }
248}
conspire/domain/fem/block/thermal/conduction/mod.rs

conspire/domain/fem/block/thermal/conduction/
mod.rs
