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 russell_lab::{Matrix, Vector};

/// Defines a line (segment) with 5 nodes (quartic functions)
///
/// ![lin5](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/test_draw_shape_lin5.svg)
///
/// # Local IDs of nodes
///
/// ```text
///   0-----3-----2-----4-----1  --> r
/// -1.0  -0.5   0.0   0.5   1.0
/// ```
///
/// # Notes
///
/// * The reference coordinates range from -1 to +1 with the geometry centred @ 0
/// * This shape is a higher-order version of [super::Lin3]
pub struct Lin5 {}

impl Lin5 {
    pub const GEO_NDIM: usize = 1;
    pub const NNODE: usize = 5;
    pub const NEDGE: usize = 0;
    pub const NFACE: usize = 0;
    pub const EDGE_NNODE: usize = 0;
    pub const FACE_NNODE: usize = 0;
    pub const FACE_NEDGE: usize = 0;
    pub const N_INTERIOR_NODE: usize = 3;

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

    pub const INTERIOR_NODES: [usize; Lin5::N_INTERIOR_NODE] = [2, 3, 4];

    /// 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 = ksi[0];

        interp[0] = (r - 1.0) * (1.0 - 2.0 * r) * r * (-1.0 - 2.0 * r) / 6.0;
        interp[1] = (1.0 - 2.0 * r) * r * (-1.0 - 2.0 * r) * (1.0 + r) / 6.0;
        interp[2] = (1.0 - r) * (1.0 - 2.0 * r) * (-1.0 - 2.0 * r) * (-1.0 - r);
        interp[3] = 4.0 * (1.0 - r) * (1.0 - 2.0 * r) * r * (-1.0 - r) / 3.0;
        interp[4] = 4.0 * (1.0 - r) * r * (-1.0 - 2.0 * r) * (-1.0 - r) / 3.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 = ksi[0];

        deriv.set(
            0,
            0,
            -((1.0 - 2.0 * r) * (r - 1.0) * r) / 3.0 - ((-2.0 * r - 1.0) * (r - 1.0) * r) / 3.0
                + ((-2.0 * r - 1.0) * (1.0 - 2.0 * r) * r) / 6.0
                + ((-2.0 * r - 1.0) * (1.0 - 2.0 * r) * (r - 1.0)) / 6.0,
        );

        deriv.set(
            1,
            0,
            -((1.0 - 2.0 * r) * r * (r + 1.0)) / 3.0 - ((-2.0 * r - 1.0) * r * (r + 1.0)) / 3.0
                + ((-2.0 * r - 1.0) * (1.0 - 2.0 * r) * (r + 1.0)) / 6.0
                + ((-2.0 * r - 1.0) * (1.0 - 2.0 * r) * r) / 6.0,
        );

        deriv.set(
            2,
            0,
            -2.0 * (1.0 - 2.0 * r) * (-r - 1.0) * (1.0 - r)
                - 2.0 * (-2.0 * r - 1.0) * (-r - 1.0) * (1.0 - r)
                - (-2.0 * r - 1.0) * (1.0 - 2.0 * r) * (1.0 - r)
                - (-2.0 * r - 1.0) * (1.0 - 2.0 * r) * (-r - 1.0),
        );

        deriv.set(
            3,
            0,
            -(8.0 * (-r - 1.0) * (1.0 - r) * r) / 3.0
                - (4.0 * (1.0 - 2.0 * r) * (1.0 - r) * r) / 3.0
                - (4.0 * (1.0 - 2.0 * r) * (-r - 1.0) * r) / 3.0
                + (4.0 * (1.0 - 2.0 * r) * (-r - 1.0) * (1.0 - r)) / 3.0,
        );

        deriv.set(
            4,
            0,
            -(8.0 * (-r - 1.0) * (1.0 - r) * r) / 3.0
                - (4.0 * (-2.0 * r - 1.0) * (1.0 - r) * r) / 3.0
                - (4.0 * (-2.0 * r - 1.0) * (-r - 1.0) * r) / 3.0
                + (4.0 * (-2.0 * r - 1.0) * (-r - 1.0) * (1.0 - r)) / 3.0,
        );
    }
}