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
use russell_lab::{Matrix, Vector};
/// Defines a line (segment) with 4 nodes (cubic functions)
///
/// 
///
/// # Local IDs of nodes
///
/// ```text
/// 0-------2-------3-------1 --> r
/// -1.0 -1/3 1/3 1.0
/// ```
///
/// # Notes
///
/// * The reference coordinates range from -1 to +1 with the geometry centred @ 0
pub struct Lin4 {}
impl Lin4 {
pub const GEO_NDIM: usize = 1;
pub const NNODE: usize = 4;
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 = 2;
#[rustfmt::skip]
pub const NODE_REFERENCE_COORDS: [[f64; Lin4::GEO_NDIM]; Lin4::NNODE] = [
[-1.0],
[ 1.0],
[-1.0 / 3.0],
[ 1.0 / 3.0],
];
pub const INTERIOR_NODES: [usize; Lin4::N_INTERIOR_NODE] = [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 = ksi[0];
interp[0] = (-9.0 * r * r * r + 9.0 * r * r + r - 1.0) / 16.0;
interp[1] = (9.0 * r * r * r + 9.0 * r * r - r - 1.0) / 16.0;
interp[2] = (27.0 * r * r * r - 9.0 * r * r - 27.0 * r + 9.0) / 16.0;
interp[3] = (-27.0 * r * r * r - 9.0 * r * r + 27.0 * r + 9.0) / 16.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 / 16.0 * (-27.0 * r * r + 18.0 * r + 1.0));
deriv.set(1, 0, 1.0 / 16.0 * (27.0 * r * r + 18.0 * r - 1.0));
deriv.set(2, 0, 1.0 / 16.0 * (81.0 * r * r - 18.0 * r - 27.0));
deriv.set(3, 0, 1.0 / 16.0 * (-81.0 * r * r - 18.0 * r + 27.0));
}
}