outram_foam_basic_lib/primitives/
vector.rs

1// SPDX-License-Identifier: GPL-3.0-only
2// Copyright (C) 2026 OUTRAM PARK contributors
3// Derived from OpenFOAM (www.openfoam.com)
4// Copyright (C) 2004-2023 OpenFOAM Foundation
5// Copyright (C) 2016-2023 OpenCFD Ltd.
6//
7// This file is part of OUTRAM PARK.
8//
9// OUTRAM PARK is free software: you can redistribute it and/or modify it
10// under the terms of the GNU General Public License as published by the
11// Free Software Foundation, either version 3 of the License, or (at your
12// option) any later version.
13//
14// OUTRAM PARK is distributed in the hope that it will be useful, but
15// WITHOUT ANY WARRANTY; without even the implied warranty of
16// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17// General Public License for more details.
18//
19// You should have received a copy of the GNU General Public License along
20// with OUTRAM PARK.  If not, see <https://www.gnu.org/licenses/>.
21
22//! 3-component vector (`Vector3`) and its OpenFOAM-style operators.
23//!
24//! Components are dimensionless `f64`; the same type is reused for positions,
25//! velocities, forces, etc. with the physical dimension carried by the caller.
26
27use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
28
29/// 3-component vector. Maps to `Foam::vector` (`Foam::Vector<scalar>`).
30/// Component layout: x, y, z.
31#[derive(Debug, Clone, Copy, PartialEq, Default)]
32pub struct Vector3 {
33    pub x: f64,
34    pub y: f64,
35    pub z: f64,
36}
37
38impl Vector3 {
39    /// Zero vector (0, 0, 0).
40    pub const ZERO: Self = Self {
41        x: 0.0,
42        y: 0.0,
43        z: 0.0,
44    };
45    /// Vector of ones (1, 1, 1).
46    pub const ONE: Self = Self {
47        x: 1.0,
48        y: 1.0,
49        z: 1.0,
50    };
51    /// Unit vector along the x axis (1, 0, 0).
52    pub const X: Self = Self {
53        x: 1.0,
54        y: 0.0,
55        z: 0.0,
56    };
57    /// Unit vector along the y axis (0, 1, 0).
58    pub const Y: Self = Self {
59        x: 0.0,
60        y: 1.0,
61        z: 0.0,
62    };
63    /// Unit vector along the z axis (0, 0, 1).
64    pub const Z: Self = Self {
65        x: 0.0,
66        y: 0.0,
67        z: 1.0,
68    };
69
70    /// Construct a vector from its x, y, z components.
71    #[inline]
72    pub fn new(x: f64, y: f64, z: f64) -> Self {
73        Self { x, y, z }
74    }
75
76    /// Squared magnitude: |v|² = x² + y² + z²
77    #[inline]
78    pub fn mag_sqr(self) -> f64 {
79        self.x * self.x + self.y * self.y + self.z * self.z
80    }
81
82    /// Magnitude: |v|
83    #[inline]
84    pub fn mag(self) -> f64 {
85        self.mag_sqr().sqrt()
86    }
87
88    /// Squared distance to another vector
89    #[inline]
90    pub fn dist_sqr(self, other: Self) -> f64 {
91        (other - self).mag_sqr()
92    }
93
94    /// Distance to another vector
95    #[inline]
96    pub fn dist(self, other: Self) -> f64 {
97        self.dist_sqr(other).sqrt()
98    }
99
100    /// Dot (inner) product. C++ `operator&(Vector, Vector)`.
101    #[inline]
102    pub fn dot(self, other: Self) -> f64 {
103        self.x * other.x + self.y * other.y + self.z * other.z
104    }
105
106    /// Cross product. C++ `operator^(Vector, Vector)`.
107    #[inline]
108    pub fn cross(self, other: Self) -> Self {
109        Self {
110            x: self.y * other.z - self.z * other.y,
111            y: self.z * other.x - self.x * other.z,
112            z: self.x * other.y - self.y * other.x,
113        }
114    }
115
116    /// Normalise to unit vector; returns zero if `|v| < tol`.
117    pub fn normalise(self, tol: f64) -> Self {
118        let s = self.mag();
119        if s < tol {
120            Self::ZERO
121        } else {
122            self / s
123        }
124    }
125
126    /// Remove the component collinear with `unit_vec`: `self - (self·unit) * unit`.
127    #[inline]
128    pub fn remove_collinear(self, unit_vec: Self) -> Self {
129        self - self.dot(unit_vec) * unit_vec
130    }
131
132    /// Linear interpolation: `(1-t)*a + t*b`.
133    #[inline]
134    pub fn lerp(a: Self, b: Self, t: f64) -> Self {
135        let ot = 1.0 - t;
136        Self {
137            x: ot * a.x + t * b.x,
138            y: ot * a.y + t * b.y,
139            z: ot * a.z + t * b.z,
140        }
141    }
142}
143
144// --- Standard arithmetic operators ---
145
146impl Neg for Vector3 {
147    type Output = Self;
148    #[inline]
149    fn neg(self) -> Self {
150        Self {
151            x: -self.x,
152            y: -self.y,
153            z: -self.z,
154        }
155    }
156}
157
158impl Add for Vector3 {
159    type Output = Self;
160    #[inline]
161    fn add(self, r: Self) -> Self {
162        Self {
163            x: self.x + r.x,
164            y: self.y + r.y,
165            z: self.z + r.z,
166        }
167    }
168}
169
170impl Sub for Vector3 {
171    type Output = Self;
172    #[inline]
173    fn sub(self, r: Self) -> Self {
174        Self {
175            x: self.x - r.x,
176            y: self.y - r.y,
177            z: self.z - r.z,
178        }
179    }
180}
181
182impl Mul<f64> for Vector3 {
183    type Output = Self;
184    #[inline]
185    fn mul(self, s: f64) -> Self {
186        Self {
187            x: self.x * s,
188            y: self.y * s,
189            z: self.z * s,
190        }
191    }
192}
193
194impl Mul<Vector3> for f64 {
195    type Output = Vector3;
196    #[inline]
197    fn mul(self, v: Vector3) -> Vector3 {
198        v * self
199    }
200}
201
202impl Div<f64> for Vector3 {
203    type Output = Self;
204    #[inline]
205    fn div(self, s: f64) -> Self {
206        Self {
207            x: self.x / s,
208            y: self.y / s,
209            z: self.z / s,
210        }
211    }
212}
213
214impl AddAssign for Vector3 {
215    #[inline]
216    fn add_assign(&mut self, r: Self) {
217        *self = *self + r;
218    }
219}
220
221impl SubAssign for Vector3 {
222    #[inline]
223    fn sub_assign(&mut self, r: Self) {
224        *self = *self - r;
225    }
226}
227
228impl MulAssign<f64> for Vector3 {
229    #[inline]
230    fn mul_assign(&mut self, s: f64) {
231        *self = *self * s;
232    }
233}
234
235impl DivAssign<f64> for Vector3 {
236    #[inline]
237    fn div_assign(&mut self, s: f64) {
238        *self = *self / s;
239    }
240}
241
242// --- Free functions mirroring OpenFOAM globals ---
243
244/// Squared magnitude |v|² of a vector.
245#[inline]
246pub fn mag_sqr(v: Vector3) -> f64 {
247    v.mag_sqr()
248}
249
250/// Magnitude |v| of a vector.
251#[inline]
252pub fn mag(v: Vector3) -> f64 {
253    v.mag()
254}
255
256/// Dot product. C++ `operator&`.
257#[inline]
258pub fn dot(a: Vector3, b: Vector3) -> f64 {
259    a.dot(b)
260}
261
262/// Cross product. C++ `operator^`.
263#[inline]
264pub fn cross(a: Vector3, b: Vector3) -> Vector3 {
265    a.cross(b)
266}
267
268/// Linear interpolation `(1-t)*a + t*b` between two vectors.
269#[inline]
270pub fn lerp(a: Vector3, b: Vector3, t: f64) -> Vector3 {
271    Vector3::lerp(a, b, t)
272}
273
274#[cfg(test)]
275mod tests {
276    use super::*;
277
278    #[test]
279    fn dot_and_cross() {
280        let x = Vector3::X;
281        let y = Vector3::Y;
282        let z = Vector3::Z;
283        assert_eq!(dot(x, x), 1.0);
284        assert_eq!(dot(x, y), 0.0);
285        assert_eq!(cross(x, y), z);
286        assert_eq!(cross(y, z), x);
287        assert_eq!(cross(z, x), y);
288    }
289
290    #[test]
291    fn cross_orthogonal_to_both_inputs() {
292        // For arbitrary vectors a, b: cross(a,b)·a == 0 and cross(a,b)·b == 0
293        let a = Vector3::new(1.0, 2.0, 3.0);
294        let b = Vector3::new(4.0, 5.0, 6.0);
295        let c = cross(a, b);
296        assert!(dot(c, a).abs() < 1e-14, "cross(a,b)·a = {}", dot(c, a));
297        assert!(dot(c, b).abs() < 1e-14, "cross(a,b)·b = {}", dot(c, b));
298    }
299
300    #[test]
301    fn mag_and_normalise() {
302        let v = Vector3::new(3.0, 4.0, 0.0);
303        assert_eq!(v.mag(), 5.0);
304        let n = v.normalise(1e-15);
305        assert!((n.mag() - 1.0).abs() < 1e-14);
306    }
307
308    #[test]
309    fn lerp() {
310        let a = Vector3::new(0.0, 0.0, 0.0);
311        let b = Vector3::new(2.0, 4.0, 6.0);
312        let m = Vector3::lerp(a, b, 0.5);
313        assert_eq!(m, Vector3::new(1.0, 2.0, 3.0));
314    }
315
316    #[test]
317    fn remove_collinear() {
318        let v = Vector3::new(3.0, 4.0, 0.0);
319        let n = v.normalise(1e-15);
320        let r = v.remove_collinear(n);
321        assert!(r.mag() < 1e-14);
322    }
323}
outram_foam_basic_lib/primitives/vector.rs

outram_foam_basic_lib/primitives/
vector.rs
