1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
macro_rules! test_thermal {
($element: ident) => {
mod thermal_element {
use super::{N, element};
use crate::{
EPSILON,
constitutive::thermal::conduction::Fourier,
fem::block::element::{
FiniteElementError,
thermal::{
ElementNodalTemperatures,
conduction::{
ElementNodalForcesThermal, ElementNodalStiffnessesThermal,
ThermalConductionFiniteElement,
},
},
},
math::assert::AssertionError,
};
mod finite_difference {
use super::*;
use $crate::units::PowerPerLengthTemperature;
const MODEL: Fourier = Fourier {
thermal_conductivity: PowerPerLengthTemperature::watts_per_meter_kelvin(1.0),
};
#[test]
fn potential() -> Result<(), AssertionError> {
let constitutive_model = MODEL;
let element = element();
let temperature = ElementNodalTemperatures::from(
[0.62895714, 0.73331084, 0.3058115, 0.08179408]
.map($crate::math::Quantity::new),
);
let mut finite_difference = $crate::math::Quantity::default();
let nodal_forces_fd: ElementNodalForcesThermal<N> = (0..N)
.map(|node| {
let mut nodal_temperatures = temperature.clone();
nodal_temperatures[node] +=
$crate::math::assert::perturbation(0.5 * EPSILON);
finite_difference =
element.potential(&constitutive_model, &nodal_temperatures)?;
nodal_temperatures[node] -= $crate::math::assert::perturbation(EPSILON);
finite_difference -=
element.potential(&constitutive_model, &nodal_temperatures)?;
// A potential per unit temperature is a power.
Ok(finite_difference
/ $crate::math::Quantity::<$crate::units::Temperature>::new(
EPSILON,
))
})
.collect::<Result<_, FiniteElementError>>()?;
$crate::math::assert::Assert::default().eq_within_fd_tol(
&nodal_forces_fd,
&element.nodal_forces(
&constitutive_model,
&ElementNodalTemperatures::from(temperature),
)?,
)
}
#[test]
fn nodal_forces() -> Result<(), AssertionError> {
let constitutive_model = MODEL;
let element = element();
let temperature = ElementNodalTemperatures::from(
[0.62895714, 0.73331084, 0.3058115, 0.08179408]
.map($crate::math::Quantity::new),
);
let mut finite_difference = $crate::math::Quantity::default();
let nodal_stiffnesses_fd: ElementNodalStiffnessesThermal<N> = (0..N)
.map(|node_a| {
(0..N)
.map(|node_b| {
let mut nodal_temperatures = temperature.clone();
nodal_temperatures[node_b] +=
$crate::math::assert::perturbation(0.5 * EPSILON);
finite_difference = element
.nodal_forces(&constitutive_model, &nodal_temperatures)?
[node_a];
nodal_temperatures[node_b] -=
$crate::math::assert::perturbation(EPSILON);
finite_difference -= element
.nodal_forces(&constitutive_model, &nodal_temperatures)?
[node_a];
// A power per unit temperature is a thermal stiffness.
Ok(finite_difference
/ $crate::math::Quantity::<$crate::units::Temperature>::new(
EPSILON,
))
})
.collect()
})
.collect::<Result<_, FiniteElementError>>()?;
$crate::math::assert::Assert::default().eq_within_fd_tol(
&nodal_stiffnesses_fd,
&element.nodal_stiffnesses(
&constitutive_model,
&ElementNodalTemperatures::from(temperature),
)?,
)
}
}
}
};
}
pub(crate) use test_thermal;