gemlab 2.0.0

Geometry and meshes laboratory for finite element analyses
Documentation 
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

use crate::integ::Gauss;
use crate::shapes::Scratchpad;
use crate::StrError;
use russell_lab::Vector;

/// Calculates the x-y-z coordinates of all integration points
///
/// This function applies the isoparametric formula
/// to all p-th integration points `ιᵖ` according to:
///
/// ```text
/// → →          → →   →
/// x(ιᵖ) = Σ Nᵐ(ξ=ιᵖ) xᵐ
///         m         
/// ```
///
/// # Input
///
/// * `integ_points` -- Integration points' constants (ngauss)
///
/// # Output
///
/// * Returns an array with `ngauss` (number of integration points) vectors, where
///   each vector has a dimension equal to `space_ndim`.
///
/// # Examples
///
/// ```
/// use gemlab::integ::Gauss;
/// use gemlab::recovery::get_points_coords;
/// use gemlab::shapes::{GeoKind, Scratchpad};
/// use gemlab::StrError;
///
/// fn main() -> Result<(), StrError> {
///     //  6 2
///     //  5 | `.    * indicates the
///     //  4 | * `.    location of ips
///     //  3 |     `.
///     //  2 |       `.
///     //  1 | *     * `.
///     //  0 0-----------1
///     //    0 1 2 3 4 5 6
///
///     let space_ndim = 2;
///     let mut pad = Scratchpad::new(space_ndim, GeoKind::Tri3)?;
///     pad.set_xx(0, 0, 0.0);
///     pad.set_xx(0, 1, 0.0);
///     pad.set_xx(1, 0, 6.0);
///     pad.set_xx(1, 1, 0.0);
///     pad.set_xx(2, 0, 0.0);
///     pad.set_xx(2, 1, 6.0);
///
///     let gauss = Gauss::new_sized(pad.kind.class(), 3)?;
///     let x_ips = get_points_coords(&mut pad, &gauss)?;
///     assert_eq!(x_ips[0].as_data(), &[1.0, 1.0]);
///     assert_eq!(x_ips[1].as_data(), &[4.0, 1.0]);
///     assert_eq!(x_ips[2].as_data(), &[1.0, 4.0]);
///     Ok(())
/// }
/// ```
pub fn get_points_coords(pad: &mut Scratchpad, gauss: &Gauss) -> Result<Vec<Vector>, StrError> {
    let space_ndim = pad.xxt.dims().0;
    let mut all_coords = Vec::new();
    let ngauss = gauss.npoint();
    for p in 0..ngauss {
        let mut x = Vector::new(space_ndim);
        pad.calc_coords(&mut x, gauss.coords(p))?;
        all_coords.push(x);
    }
    Ok(all_coords)
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod tests {
    use super::get_points_coords;
    use crate::integ::Gauss;
    use crate::shapes::{GeoKind, Scratchpad};
    use russell_lab::approx_eq;

    #[test]
    pub fn points_coords_works() {
        //  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 (w, h) = (20.0, 10.0);
        let space_ndim = 2;
        let mut pad = Scratchpad::new(space_ndim, GeoKind::Qua4).unwrap();
        pad.set_xx(0, 0, 0.0);
        pad.set_xx(0, 1, 0.0);
        pad.set_xx(1, 0, w);
        pad.set_xx(1, 1, 0.0);
        pad.set_xx(2, 0, w);
        pad.set_xx(2, 1, h);
        pad.set_xx(3, 0, 0.0);
        pad.set_xx(3, 1, h);

        let gauss = Gauss::new_sized(pad.kind.class(), 4).unwrap();
        let x_ips = get_points_coords(&mut pad, &gauss).unwrap();

        approx_eq(x_ips[0][0], w * (1.0 - f64::sqrt(3.0) / 3.0) / 2.0, 1e-15);
        approx_eq(x_ips[0][1], h * (1.0 - f64::sqrt(3.0) / 3.0) / 2.0, 1e-15);
        approx_eq(x_ips[1][0], w * (1.0 + f64::sqrt(3.0) / 3.0) / 2.0, 1e-15);
        approx_eq(x_ips[1][1], h * (1.0 - f64::sqrt(3.0) / 3.0) / 2.0, 1e-15);
        approx_eq(x_ips[2][0], w * (1.0 - f64::sqrt(3.0) / 3.0) / 2.0, 1e-15);
        approx_eq(x_ips[2][1], h * (1.0 + f64::sqrt(3.0) / 3.0) / 2.0, 1e-15);
        approx_eq(x_ips[3][0], w * (1.0 + f64::sqrt(3.0) / 3.0) / 2.0, 1e-15);
        approx_eq(x_ips[3][1], h * (1.0 + f64::sqrt(3.0) / 3.0) / 2.0, 1e-15);
    }
}