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}