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