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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
use super::Scratchpad;
use crate::StrError;
use russell_lab::{mat_vec_mul, Vector};
impl Scratchpad {
/// Calculates the real coordinates x from reference coordinates ξ
///
/// This function uses the isoparametric formula to calculate x given ξ:
///
/// ```text
/// → → → →
/// x(ξ) = Σ Nᵐ(ξ) xᵐ
/// m
///
/// x := Xᵀ ⋅ N
/// ```
///
/// # Output
///
/// * `x` -- real coordinates (space_ndim)
/// * `pad.interp` -- (nnode) interpolation functions @ ξ
///
/// # Input
///
/// * `ksi` -- reference coordinates ξ with len ≥ geo_ndim
///
/// # Examples
///
/// ```
/// use gemlab::shapes::{GeoKind, Scratchpad};
/// use gemlab::StrError;
/// use russell_lab::{Vector, vec_approx_eq};
///
/// fn main() -> Result<(), StrError> {
/// // 3-------------2 ξ₀ ξ₁
/// // | ξ₁ | node r s
/// // | | | 0 -1.0 -1.0
/// // | +--ξ₀ | 1 1.0 -1.0
/// // | | 2 1.0 1.0
/// // | | 3 -1.0 1.0
/// // 0-------------1
///
/// let (x0, y0) = (3.0, 4.0);
/// let (w, h) = (10.0, 5.0);
/// let space_ndim = 2;
/// let mut pad = Scratchpad::new(space_ndim, GeoKind::Qua4)?;
/// pad.set_xx(0, 0, x0);
/// pad.set_xx(0, 1, y0);
/// pad.set_xx(1, 0, x0 + w);
/// pad.set_xx(1, 1, y0);
/// pad.set_xx(2, 0, x0 + w);
/// pad.set_xx(2, 1, y0 + h);
/// pad.set_xx(3, 0, x0);
/// pad.set_xx(3, 1, y0 + h);
///
/// let mut x = Vector::new(2);
/// pad.calc_coords(&mut x, &[0.0, 0.0])?;
/// vec_approx_eq(&x, &[x0 + w / 2.0, y0 + h / 2.0], 1e-15);
/// Ok(())
/// }
/// ```
pub fn calc_coords(&mut self, x: &mut Vector, ksi: &[f64]) -> Result<(), StrError> {
if !self.ok_xxt {
return Err("all components of the coordinates matrix must be set first");
}
let space_ndim = self.jacobian.dims().0;
if x.dim() != space_ndim {
return Err("x.dim() must be equal to space_ndim");
}
self.calc_interp(ksi);
mat_vec_mul(x, 1.0, &self.xxt, &self.interp).unwrap();
Ok(())
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
#[cfg(test)]
mod tests {
use crate::shapes::scratchpad_testing::aux;
use crate::shapes::{GeoKind, Scratchpad};
use russell_lab::math::ONE_BY_3;
use russell_lab::{vec_approx_eq, Vector};
#[test]
fn calc_coords_handles_errors() {
let mut x = Vector::new(1);
let mut pad = Scratchpad::new(2, GeoKind::Tri3).unwrap();
assert_eq!(
pad.calc_coords(&mut x, &[0.0, 0.0]).err(),
Some("all components of the coordinates matrix must be set first")
);
pad.set_xx(2, 1, 0.0); // setting the last component
// (cannot really check that all components have been set)
assert_eq!(
pad.calc_coords(&mut x, &[0.0, 0.0]).err(),
Some("x.dim() must be equal to space_ndim")
);
}
#[test]
fn calc_coords_works() {
// kind, tol, tol_in for the case inside the shape
let problem = vec![
(GeoKind::Tri3, 1e-15, 0.35), // linear maps are inaccurate for the circular wedge
(GeoKind::Tri6, 1e-15, 0.013), // << quadratic mapping is inaccurate as well
(GeoKind::Tri10, 1e-14, 1e-14),
(GeoKind::Tri15, 1e-14, 1e-5), // << this triangle is inaccurate as well here
(GeoKind::Qua4, 1e-15, 0.19), // linear maps are inaccurate for the circular wedge
(GeoKind::Qua8, 1e-15, 1e-14),
(GeoKind::Qua17, 1e-15, 1e-15),
(GeoKind::Tet4, 1e-15, 0.35), // linear tetrahedron is also inaccurate here
(GeoKind::Tet10, 1e-15, 0.013), // quadratic tetrahedron is also inaccurate here
(GeoKind::Tet20, 1e-14, 1e-14), // cubic tetrahedron
(GeoKind::Hex8, 1e-14, 0.19), // bi-linear maps are inaccurate for the circular wedge
(GeoKind::Hex20, 1e-15, 1e-15),
(GeoKind::Hex32, 1e-15, 0.00012), // TODO: check why this tolerance is high
];
// loop over shapes
for (kind, tol, tol_in) in problem {
println!("kind = {:?}", kind);
// scratchpad with coordinates
let geo_ndim = kind.ndim();
let space_ndim = usize::max(2, geo_ndim);
let mut pad = aux::gen_scratchpad_with_coords(space_ndim, kind);
// loop over nodes of shape
let nnode = kind.nnode();
let mut x = Vector::new(space_ndim);
let mut x_correct = Vector::new(space_ndim);
let (ksi_min, ksi_del) = kind.ksi_min_ksi_del();
for m in 0..nnode {
// get ξᵐ corresponding to node m
let ksi = kind.reference_coords(m);
// calculate xᵐ(ξᵐ) using the isoparametric formula
pad.calc_coords(&mut x, ksi).unwrap();
// compare xᵐ with generated coordinates
aux::map_point_coords(&mut x_correct, ksi, ksi_min, ksi_del);
vec_approx_eq(&x, &x_correct, tol);
}
// test again inside the reference domain
let ksi_in = if kind.is_tri_or_tet() {
vec![ONE_BY_3; geo_ndim]
} else {
vec![0.0; geo_ndim]
};
pad.calc_coords(&mut x, &ksi_in).unwrap();
aux::map_point_coords(&mut x_correct, &ksi_in, ksi_min, ksi_del);
vec_approx_eq(&x, &x_correct, tol_in);
}
}
}