outram_foam_basic_lib/primitives/
spherical_tensor.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
22use std::ops::{Add, Div, Mul, Neg, Sub};
23
24/// Isotropic diagonal tensor: represents `ii * I` where `I` is the 3×3 identity.
25/// Maps to `Foam::SphericalTensor<scalar>` (`SphericalTensorI.H`).
26#[derive(Debug, Clone, Copy, PartialEq, Default)]
27pub struct SphericalTensor {
28    pub ii: f64,
29}
30
31impl SphericalTensor {
32    pub const ZERO: Self = Self { ii: 0.0 };
33    pub const IDENTITY: Self = Self { ii: 1.0 };
34
35    #[inline]
36    pub fn new(ii: f64) -> Self {
37        Self { ii }
38    }
39
40    /// Trace = 3 * ii
41    #[inline]
42    pub fn tr(self) -> f64 {
43        3.0 * self.ii
44    }
45
46    /// Frobenius norm squared = 3 * ii²
47    #[inline]
48    pub fn mag_sqr(self) -> f64 {
49        3.0 * self.ii * self.ii
50    }
51
52    #[inline]
53    pub fn mag(self) -> f64 {
54        self.mag_sqr().sqrt()
55    }
56
57    /// Diagonal norm squared (sum of squared diagonal entries = 3*ii²)
58    #[inline]
59    pub fn diag_sqr(self) -> f64 {
60        3.0 * self.ii * self.ii
61    }
62
63    /// Determinant = ii³
64    #[inline]
65    pub fn det(self) -> f64 {
66        self.ii * self.ii * self.ii
67    }
68
69    /// Inverse: SphericalTensor(1/ii)
70    #[inline]
71    pub fn inv(self) -> Self {
72        Self { ii: 1.0 / self.ii }
73    }
74
75    /// Double inner-product with itself: 3 * ii²
76    #[inline]
77    pub fn double_inner(self, rhs: Self) -> f64 {
78        3.0 * self.ii * rhs.ii
79    }
80
81    /// Linear interpolation
82    #[inline]
83    pub fn lerp(a: Self, b: Self, t: f64) -> Self {
84        Self {
85            ii: (1.0 - t) * a.ii + t * b.ii,
86        }
87    }
88}
89
90// --- Arithmetic operators ---
91
92impl Neg for SphericalTensor {
93    type Output = Self;
94    #[inline]
95    fn neg(self) -> Self {
96        Self { ii: -self.ii }
97    }
98}
99
100impl Add for SphericalTensor {
101    type Output = Self;
102    #[inline]
103    fn add(self, rhs: Self) -> Self {
104        Self {
105            ii: self.ii + rhs.ii,
106        }
107    }
108}
109
110impl Sub for SphericalTensor {
111    type Output = Self;
112    #[inline]
113    fn sub(self, rhs: Self) -> Self {
114        Self {
115            ii: self.ii - rhs.ii,
116        }
117    }
118}
119
120impl Mul<f64> for SphericalTensor {
121    type Output = Self;
122    #[inline]
123    fn mul(self, s: f64) -> Self {
124        Self { ii: self.ii * s }
125    }
126}
127
128impl Mul<SphericalTensor> for f64 {
129    type Output = SphericalTensor;
130    #[inline]
131    fn mul(self, st: SphericalTensor) -> SphericalTensor {
132        st * self
133    }
134}
135
136impl Div<f64> for SphericalTensor {
137    type Output = Self;
138    #[inline]
139    fn div(self, s: f64) -> Self {
140        Self { ii: self.ii / s }
141    }
142}
143
144/// `scalar / SphericalTensor` — maps to C++ `operator/(Cmpt, SphericalTensor)`
145impl Div<SphericalTensor> for f64 {
146    type Output = SphericalTensor;
147    #[inline]
148    fn div(self, st: SphericalTensor) -> SphericalTensor {
149        SphericalTensor { ii: self / st.ii }
150    }
151}
152
153// --- Free functions mirroring OpenFOAM globals ---
154
155#[inline]
156pub fn tr(st: SphericalTensor) -> f64 {
157    st.tr()
158}
159
160#[inline]
161pub fn det(st: SphericalTensor) -> f64 {
162    st.det()
163}
164
165#[inline]
166pub fn inv(st: SphericalTensor) -> SphericalTensor {
167    st.inv()
168}
169
170#[inline]
171pub fn mag_sqr(st: SphericalTensor) -> f64 {
172    st.mag_sqr()
173}
174
175#[inline]
176pub fn lerp(a: SphericalTensor, b: SphericalTensor, t: f64) -> SphericalTensor {
177    SphericalTensor::lerp(a, b, t)
178}
179
180#[cfg(test)]
181mod tests {
182    use super::*;
183
184    #[test]
185    fn identity_properties() {
186        let i = SphericalTensor::IDENTITY;
187        assert_eq!(i.tr(), 3.0);
188        assert_eq!(i.det(), 1.0);
189        assert_eq!(i.inv().ii, 1.0);
190        assert_eq!(i.mag_sqr(), 3.0);
191    }
192
193    #[test]
194    fn arithmetic() {
195        let a = SphericalTensor::new(2.0);
196        let b = SphericalTensor::new(3.0);
197        assert_eq!((a + b).ii, 5.0);
198        assert_eq!((b - a).ii, 1.0);
199        assert_eq!((a * 3.0).ii, 6.0);
200        assert_eq!((3.0 * a).ii, 6.0);
201        assert_eq!((a / 2.0).ii, 1.0);
202    }
203
204    #[test]
205    fn double_inner() {
206        let a = SphericalTensor::new(2.0);
207        let b = SphericalTensor::new(3.0);
208        assert_eq!(a.double_inner(b), 3.0 * 2.0 * 3.0);
209    }
210}
outram_foam_basic_lib/primitives/spherical_tensor.rs

outram_foam_basic_lib/primitives/
spherical_tensor.rs
