outram-foam-basic-lib 0.1.0

OpenFOAM primitive functions translated to Rust -- an independent OUTRAM PARK fork, not the official OpenFOAM software
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

// SPDX-License-Identifier: GPL-3.0-only
// Copyright (C) 2026 OUTRAM PARK contributors
// Derived from OpenFOAM (www.openfoam.com)
// Copyright (C) 2004-2023 OpenFOAM Foundation
// Copyright (C) 2016-2023 OpenCFD Ltd.
//
// This file is part of OUTRAM PARK.
//
// OUTRAM PARK is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the
// Free Software Foundation, either version 3 of the License, or (at your
// option) any later version.
//
// OUTRAM PARK is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with OUTRAM PARK.  If not, see <https://www.gnu.org/licenses/>.

/// Linear 1-D interpolation over a sorted table `(xs, ys)`.
///
/// Clamps to the endpoint values outside the table range.
/// Assumes `xs` is sorted in ascending order.
/// Maps to `Foam::interpolateXY(scalar, UList<scalar>&, UList<Type>&)`.
pub fn interpolate_xy(x: f64, xs: &[f64], ys: &[f64]) -> f64 {
    let n = xs.len();
    assert_eq!(n, ys.len(), "xs and ys must have the same length");

    if n == 0 {
        return 0.0;
    }
    if n == 1 || x <= xs[0] {
        return ys[0];
    }
    if x >= xs[n - 1] {
        return ys[n - 1];
    }

    // Binary search: hi = first index where xs[hi] >= x
    let hi = xs.partition_point(|&v| v < x);
    let lo = hi - 1;

    let t = (x - xs[lo]) / (xs[hi] - xs[lo]);
    ys[lo] + t * (ys[hi] - ys[lo])
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn exact_at_knots() {
        let xs = [0.0, 1.0, 2.0, 3.0];
        let ys = [0.0, 1.0, 4.0, 9.0];
        for i in 0..4 {
            let v = interpolate_xy(xs[i], &xs, &ys);
            assert!((v - ys[i]).abs() < 1e-14, "at {}: got {}", xs[i], v);
        }
    }

    #[test]
    fn midpoints() {
        let xs = [0.0, 1.0, 2.0];
        let ys = [0.0, 2.0, 6.0];
        // 0.5 → 1.0
        let v = interpolate_xy(0.5, &xs, &ys);
        assert!((v - 1.0).abs() < 1e-14, "got {}", v);
        // 1.5 → 4.0
        let v = interpolate_xy(1.5, &xs, &ys);
        assert!((v - 4.0).abs() < 1e-14, "got {}", v);
    }

    #[test]
    fn clamp_left() {
        let xs = [1.0, 2.0, 3.0];
        let ys = [10.0, 20.0, 30.0];
        assert_eq!(interpolate_xy(-5.0, &xs, &ys), 10.0);
        assert_eq!(interpolate_xy(1.0, &xs, &ys), 10.0);
    }

    #[test]
    fn clamp_right() {
        let xs = [1.0, 2.0, 3.0];
        let ys = [10.0, 20.0, 30.0];
        assert_eq!(interpolate_xy(3.0, &xs, &ys), 30.0);
        assert_eq!(interpolate_xy(100.0, &xs, &ys), 30.0);
    }

    #[test]
    fn single_point() {
        let xs = [5.0];
        let ys = [42.0];
        assert_eq!(interpolate_xy(0.0, &xs, &ys), 42.0);
        assert_eq!(interpolate_xy(5.0, &xs, &ys), 42.0);
        assert_eq!(interpolate_xy(9.0, &xs, &ys), 42.0);
    }

    #[test]
    fn two_points() {
        let xs = [0.0, 2.0];
        let ys = [0.0, 4.0];
        let v = interpolate_xy(1.0, &xs, &ys);
        assert!((v - 2.0).abs() < 1e-14);
    }
}