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
117
118
119
120
121
122
123

use russell_lab::{Matrix, Vector};

/// Defines a quadrilateral with 4 nodes (linear edges)
///
/// ![qua4](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/test_draw_shape_qua4.svg)
///
/// # Local IDs of nodes
///
/// ```text
/// 3-----------2
/// |     s     |              r     s
/// |     |     |      p:0 [-1.0, -1.0]
/// |     +--r  |      p:1 [ 1.0, -1.0]
/// |           |      p:2 [ 1.0,  1.0]
/// |           |      p:3 [-1.0,  1.0]
/// 0-----------1
/// ```
///
/// # Local IDs of edges
///
/// ```text
///        2
///  3-----------2         p0  p1
///  |           |     e:0 [1, 0]
///  |           |     e:1 [2, 1]
/// 3|           |1    e:2 [3, 2]
///  |           |     e:3 [0, 3]
///  |           |
///  0-----------1
///        0
/// ```
///
/// # Notes
///
/// * The reference coordinates range from -1 to +1 with the geometry centred @ 0
/// * The order of edge nodes is such that the normals are outward
/// * The order of edge nodes corresponds to [super::Lin2] nodes
pub struct Qua4 {}

impl Qua4 {
    pub const GEO_NDIM: usize = 2;
    pub const NNODE: usize = 4;
    pub const NEDGE: usize = 4;
    pub const NFACE: usize = 0;
    pub const EDGE_NNODE: usize = 2;
    pub const FACE_NNODE: usize = 0;
    pub const FACE_NEDGE: usize = 0;
    pub const TRIANGULATE_NTRIANGLE: usize = 2;

    #[rustfmt::skip]
    pub const EDGE_NODE_IDS: [[usize; Qua4::EDGE_NNODE]; Qua4::NEDGE] = [
        [1, 0],
        [2, 1],
        [3, 2],
        [0, 3],
    ];

    #[rustfmt::skip]
    pub const EDGE_NODE_IDS_INWARD: [[usize; Qua4::EDGE_NNODE]; Qua4::NEDGE] = [
        [0, 1],
        [1, 2],
        [2, 3],
        [3, 0],
    ];

    #[rustfmt::skip]
    pub const NODE_REFERENCE_COORDS: [[f64; Qua4::GEO_NDIM]; Qua4::NNODE] = [
        [-1.0, -1.0],
        [ 1.0, -1.0],
        [ 1.0,  1.0],
        [-1.0,  1.0],
    ];

    #[rustfmt::skip]
    pub const TRIANGULATE_TRIANGLES: [[usize; 3]; Qua4::TRIANGULATE_NTRIANGLE] = [
        [0, 1, 3],
        [1, 2, 3],
    ];

    /// Computes the interpolation functions
    ///
    /// # Output
    ///
    /// * `interp` -- interpolation function evaluated at ksi (nnode)
    ///
    /// # Input
    ///
    /// * `ksi` -- reference coordinates with length ≥ geo_ndim
    pub fn calc_interp(interp: &mut Vector, ksi: &[f64]) {
        let (r, s) = (ksi[0], ksi[1]);

        interp[0] = (1.0 - r - s + r * s) / 4.0;
        interp[1] = (1.0 + r - s - r * s) / 4.0;
        interp[2] = (1.0 + r + s + r * s) / 4.0;
        interp[3] = (1.0 - r + s - r * s) / 4.0;
    }

    /// Computes the derivatives of interpolation functions with respect to the reference coordinates
    ///
    /// # Output
    ///
    /// * `deriv` -- derivatives of the interpolation function with respect to
    ///   the reference coordinates ksi, evaluated at ksi (nnode,geo_ndim)
    ///
    /// # Input
    ///
    /// * `ksi` -- reference coordinates with length ≥ geo_ndim
    pub fn calc_deriv(deriv: &mut Matrix, ksi: &[f64]) {
        let (r, s) = (ksi[0], ksi[1]);

        deriv.set(0, 0, (-1.0 + s) / 4.0);
        deriv.set(0, 1, (-1.0 + r) / 4.0);

        deriv.set(1, 0, (1.0 - s) / 4.0);
        deriv.set(1, 1, (-1.0 - r) / 4.0);

        deriv.set(2, 0, (1.0 + s) / 4.0);
        deriv.set(2, 1, (1.0 + r) / 4.0);

        deriv.set(3, 0, (-1.0 - s) / 4.0);
        deriv.set(3, 1, (1.0 - r) / 4.0);
    }
}