outram_foam_basic_lib/primitives/
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
22//! Full (non-symmetric) 3×3 tensor (`Tensor`) and its OpenFOAM-style
23//! operators, invariants, and decompositions.
24//!
25//! Components are dimensionless `f64` stored row-major
26//! (xx, xy, xz, yx, yy, yz, zx, zy, zz). Operator names follow OpenFOAM:
27//! `mat_mul`/`mat_vec` are the single inner product (`&`), `double_inner`
28//! the double contraction (`&&`), and the dyadic/outer product is `Mul`
29//! (`*`) of two vectors.
30
31use std::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign};
32
33use super::spherical_tensor::SphericalTensor;
34use super::symm_tensor::SymmTensor;
35use super::vector::Vector3;
36
37/// Full (non-symmetric) 3×3 tensor stored row-major.
38/// Component order: xx, xy, xz, yx, yy, yz, zx, zy, zz.
39/// Maps to `Foam::tensor` (`Foam::Tensor<scalar>`).
40#[derive(Debug, Clone, Copy, PartialEq, Default)]
41pub struct Tensor {
42    pub xx: f64,
43    pub xy: f64,
44    pub xz: f64,
45    pub yx: f64,
46    pub yy: f64,
47    pub yz: f64,
48    pub zx: f64,
49    pub zy: f64,
50    pub zz: f64,
51}
52
53impl Tensor {
54    /// The zero tensor (all nine components 0).
55    pub const ZERO: Self = Self {
56        xx: 0.0,
57        xy: 0.0,
58        xz: 0.0,
59        yx: 0.0,
60        yy: 0.0,
61        yz: 0.0,
62        zx: 0.0,
63        zy: 0.0,
64        zz: 0.0,
65    };
66    /// The identity tensor (unit diagonal, zero off-diagonal).
67    pub const IDENTITY: Self = Self {
68        xx: 1.0,
69        xy: 0.0,
70        xz: 0.0,
71        yx: 0.0,
72        yy: 1.0,
73        yz: 0.0,
74        zx: 0.0,
75        zy: 0.0,
76        zz: 1.0,
77    };
78
79    /// Construct a tensor from its nine components in row-major order.
80    #[allow(clippy::too_many_arguments)]
81    #[inline]
82    pub fn new(
83        xx: f64,
84        xy: f64,
85        xz: f64,
86        yx: f64,
87        yy: f64,
88        yz: f64,
89        zx: f64,
90        zy: f64,
91        zz: f64,
92    ) -> Self {
93        Self {
94            xx,
95            xy,
96            xz,
97            yx,
98            yy,
99            yz,
100            zx,
101            zy,
102            zz,
103        }
104    }
105
106    /// Construct from three row vectors.
107    #[inline]
108    pub fn from_rows(x: Vector3, y: Vector3, z: Vector3) -> Self {
109        Self {
110            xx: x.x,
111            xy: x.y,
112            xz: x.z,
113            yx: y.x,
114            yy: y.y,
115            yz: y.z,
116            zx: z.x,
117            zy: z.y,
118            zz: z.z,
119        }
120    }
121
122    /// Construct from three column vectors.
123    #[inline]
124    pub fn from_cols(x: Vector3, y: Vector3, z: Vector3) -> Self {
125        Self {
126            xx: x.x,
127            xy: y.x,
128            xz: z.x,
129            yx: x.y,
130            yy: y.y,
131            yz: z.y,
132            zx: x.z,
133            zy: y.z,
134            zz: z.z,
135        }
136    }
137
138    /// First (x) row as a vector: (xx, xy, xz).
139    #[inline]
140    pub fn row_x(self) -> Vector3 {
141        Vector3::new(self.xx, self.xy, self.xz)
142    }
143    /// Second (y) row as a vector: (yx, yy, yz).
144    #[inline]
145    pub fn row_y(self) -> Vector3 {
146        Vector3::new(self.yx, self.yy, self.yz)
147    }
148    /// Third (z) row as a vector: (zx, zy, zz).
149    #[inline]
150    pub fn row_z(self) -> Vector3 {
151        Vector3::new(self.zx, self.zy, self.zz)
152    }
153
154    /// First (x) column as a vector: (xx, yx, zx).
155    #[inline]
156    pub fn col_x(self) -> Vector3 {
157        Vector3::new(self.xx, self.yx, self.zx)
158    }
159    /// Second (y) column as a vector: (xy, yy, zy).
160    #[inline]
161    pub fn col_y(self) -> Vector3 {
162        Vector3::new(self.xy, self.yy, self.zy)
163    }
164    /// Third (z) column as a vector: (xz, yz, zz).
165    #[inline]
166    pub fn col_z(self) -> Vector3 {
167        Vector3::new(self.xz, self.yz, self.zz)
168    }
169
170    /// Diagonal as a vector
171    #[inline]
172    pub fn diag(self) -> Vector3 {
173        Vector3::new(self.xx, self.yy, self.zz)
174    }
175
176    /// Trace
177    #[inline]
178    pub fn tr(self) -> f64 {
179        self.xx + self.yy + self.zz
180    }
181
182    /// Sum of squared diagonal entries (not Frobenius)
183    #[inline]
184    pub fn diag_sqr(self) -> f64 {
185        self.xx * self.xx + self.yy * self.yy + self.zz * self.zz
186    }
187
188    /// Transpose. C++ `.T()`.
189    #[inline]
190    pub fn transpose(self) -> Self {
191        Self {
192            xx: self.xx,
193            xy: self.yx,
194            xz: self.zx,
195            yx: self.xy,
196            yy: self.yy,
197            yz: self.zy,
198            zx: self.xz,
199            zy: self.yz,
200            zz: self.zz,
201        }
202    }
203
204    /// Determinant
205    #[inline]
206    pub fn det(self) -> f64 {
207        self.xx * (self.yy * self.zz - self.yz * self.zy)
208            + self.xy * (self.yz * self.zx - self.yx * self.zz)
209            + self.xz * (self.yx * self.zy - self.yy * self.zx)
210    }
211
212    /// Adjunct (transpose of cofactor matrix)
213    #[inline]
214    pub fn adjunct(self) -> Self {
215        Self {
216            xx: self.yy * self.zz - self.zy * self.yz,
217            xy: self.xz * self.zy - self.xy * self.zz,
218            xz: self.xy * self.yz - self.xz * self.yy,
219
220            yx: self.zx * self.yz - self.yx * self.zz,
221            yy: self.xx * self.zz - self.xz * self.zx,
222            yz: self.yx * self.xz - self.xx * self.yz,
223
224            zx: self.yx * self.zy - self.yy * self.zx,
225            zy: self.xy * self.zx - self.xx * self.zy,
226            zz: self.xx * self.yy - self.yx * self.xy,
227        }
228    }
229
230    /// Cofactor matrix = adjunct().T()
231    #[inline]
232    pub fn cof(self) -> Self {
233        self.adjunct().transpose()
234    }
235
236    /// Inverse = adjunct / det. Panics (debug) if singular.
237    #[inline]
238    pub fn inv(self) -> Self {
239        let d = self.det();
240        debug_assert!(d.abs() > 0.0, "Tensor is singular");
241        self.adjunct() / d
242    }
243
244    /// Inverse with 2-D fallback: returns ZERO if nearly singular.
245    pub fn safe_inv(self) -> Self {
246        use super::scalar::{ROOT_VSMALL, SMALL};
247        let diag_sqr = self.diag_sqr();
248        let threshold = SMALL * diag_sqr;
249        let mut work = self;
250        let small_xx = self.xx * self.xx < threshold;
251        let small_yy = self.yy * self.yy < threshold;
252        let small_zz = self.zz * self.zz < threshold;
253        if small_xx {
254            work.xx += 1.0;
255        }
256        if small_yy {
257            work.yy += 1.0;
258        }
259        if small_zz {
260            work.zz += 1.0;
261        }
262        let d = work.det();
263        if d.abs() < ROOT_VSMALL {
264            return Self::ZERO;
265        }
266        let mut result = work.adjunct() / d;
267        if small_xx {
268            result.xx -= 1.0;
269        }
270        if small_yy {
271            result.yy -= 1.0;
272        }
273        if small_zz {
274            result.zz -= 1.0;
275        }
276        result
277    }
278
279    /// Matrix multiply: `self & rhs`. C++ `operator&(Tensor, Tensor)` / `.inner(t2)`.
280    #[inline]
281    pub fn mat_mul(self, t: Self) -> Self {
282        Self {
283            xx: self.xx * t.xx + self.xy * t.yx + self.xz * t.zx,
284            xy: self.xx * t.xy + self.xy * t.yy + self.xz * t.zy,
285            xz: self.xx * t.xz + self.xy * t.yz + self.xz * t.zz,
286
287            yx: self.yx * t.xx + self.yy * t.yx + self.yz * t.zx,
288            yy: self.yx * t.xy + self.yy * t.yy + self.yz * t.zy,
289            yz: self.yx * t.xz + self.yy * t.yz + self.yz * t.zz,
290
291            zx: self.zx * t.xx + self.zy * t.yx + self.zz * t.zx,
292            zy: self.zx * t.xy + self.zy * t.yy + self.zz * t.zy,
293            zz: self.zx * t.xz + self.zy * t.yz + self.zz * t.zz,
294        }
295    }
296
297    /// Element-wise product (Schur/Hadamard product).
298    #[inline]
299    pub fn schur(self, t: Self) -> Self {
300        Self {
301            xx: self.xx * t.xx,
302            xy: self.xy * t.xy,
303            xz: self.xz * t.xz,
304            yx: self.yx * t.yx,
305            yy: self.yy * t.yy,
306            yz: self.yz * t.yz,
307            zx: self.zx * t.zx,
308            zy: self.zy * t.zy,
309            zz: self.zz * t.zz,
310        }
311    }
312
313    /// Matrix-vector multiply: `T · v`. C++ `operator&(Tensor, Vector)`.
314    #[inline]
315    pub fn mat_vec(self, v: Vector3) -> Vector3 {
316        Vector3::new(
317            self.xx * v.x + self.xy * v.y + self.xz * v.z,
318            self.yx * v.x + self.yy * v.y + self.yz * v.z,
319            self.zx * v.x + self.zy * v.y + self.zz * v.z,
320        )
321    }
322
323    /// Vector-matrix multiply: `v · T`. C++ `operator&(Vector, Tensor)`.
324    #[inline]
325    pub fn vec_mat(v: Vector3, t: Self) -> Vector3 {
326        Vector3::new(
327            v.x * t.xx + v.y * t.yx + v.z * t.zx,
328            v.x * t.xy + v.y * t.yy + v.z * t.zy,
329            v.x * t.xz + v.y * t.yz + v.z * t.zz,
330        )
331    }
332
333    /// Double contraction (full Frobenius inner product). C++ `operator&&(Tensor, Tensor)`.
334    #[inline]
335    pub fn double_inner(self, t: Self) -> f64 {
336        self.xx * t.xx
337            + self.xy * t.xy
338            + self.xz * t.xz
339            + self.yx * t.yx
340            + self.yy * t.yy
341            + self.yz * t.yz
342            + self.zx * t.zx
343            + self.zy * t.zy
344            + self.zz * t.zz
345    }
346
347    /// Symmetric part: `0.5*(T + T^T)`. Returns `SymmTensor`.
348    #[inline]
349    pub fn symm(self) -> SymmTensor {
350        SymmTensor {
351            xx: self.xx,
352            xy: 0.5 * (self.xy + self.yx),
353            xz: 0.5 * (self.xz + self.zx),
354            yy: self.yy,
355            yz: 0.5 * (self.yz + self.zy),
356            zz: self.zz,
357        }
358    }
359
360    /// Twice the symmetric part: `T + T^T`. Returns `SymmTensor`.
361    #[inline]
362    pub fn two_symm(self) -> SymmTensor {
363        SymmTensor {
364            xx: 2.0 * self.xx,
365            xy: self.xy + self.yx,
366            xz: self.xz + self.zx,
367            yy: 2.0 * self.yy,
368            yz: self.yz + self.zy,
369            zz: 2.0 * self.zz,
370        }
371    }
372
373    /// Skew-symmetric (antisymmetric) part: `0.5*(T - T^T)`.
374    #[inline]
375    pub fn skew(self) -> Self {
376        Self {
377            xx: 0.0,
378            xy: 0.5 * (self.xy - self.yx),
379            xz: 0.5 * (self.xz - self.zx),
380            yx: 0.5 * (self.yx - self.xy),
381            yy: 0.0,
382            yz: 0.5 * (self.yz - self.zy),
383            zx: 0.5 * (self.zx - self.xz),
384            zy: 0.5 * (self.zy - self.yz),
385            zz: 0.0,
386        }
387    }
388
389    /// Deviatoric part: `T - (tr/3)*I`.
390    #[inline]
391    pub fn dev(self) -> Self {
392        self - SphericalTensor::new(self.tr() / 3.0)
393    }
394
395    /// Two-thirds deviatoric: `T - (2*tr/3)*I`.
396    #[inline]
397    pub fn dev2(self) -> Self {
398        self - SphericalTensor::new(2.0 / 3.0 * self.tr())
399    }
400
401    /// Deviatoric of symmetric part: `symm(T) - (tr/3)*I`. Returns `SymmTensor`.
402    #[inline]
403    pub fn dev_symm(self) -> SymmTensor {
404        let sph = self.tr() / 3.0;
405        SymmTensor {
406            xx: self.xx - sph,
407            xy: 0.5 * (self.xy + self.yx),
408            xz: 0.5 * (self.xz + self.zx),
409            yy: self.yy - sph,
410            yz: 0.5 * (self.yz + self.zy),
411            zz: self.zz - sph,
412        }
413    }
414
415    /// Deviatoric of twice the symmetric part: `twoSymm(T) - (2*tr/3)*I`. Returns `SymmTensor`.
416    #[inline]
417    pub fn dev_two_symm(self) -> SymmTensor {
418        let sph = 2.0 / 3.0 * self.tr();
419        SymmTensor {
420            xx: 2.0 * self.xx - sph,
421            xy: self.xy + self.yx,
422            xz: self.xz + self.zx,
423            yy: 2.0 * self.yy - sph,
424            yz: self.yz + self.zy,
425            zz: 2.0 * self.zz - sph,
426        }
427    }
428
429    /// Hodge dual as a Vector. C++ `operator*(Tensor)`.
430    #[inline]
431    pub fn hodge_dual(self) -> Vector3 {
432        Vector3::new(self.yz, -self.xz, self.xy)
433    }
434
435    /// First invariant: trace
436    #[inline]
437    pub fn invariant_i(self) -> f64 {
438        self.tr()
439    }
440
441    /// Second invariant: (xx*yy + yy*zz + xx*zz) - (xy*yx + yz*zy + xz*zx)
442    #[inline]
443    pub fn invariant_ii(self) -> f64 {
444        self.xx * self.yy + self.yy * self.zz + self.xx * self.zz
445            - self.xy * self.yx
446            - self.yz * self.zy
447            - self.xz * self.zx
448    }
449
450    /// Third invariant: determinant
451    #[inline]
452    pub fn invariant_iii(self) -> f64 {
453        self.det()
454    }
455
456    /// True if approximately the identity.
457    pub fn is_identity(self, tol: f64) -> bool {
458        (self.xx - 1.0).abs() < tol
459            && (self.yy - 1.0).abs() < tol
460            && (self.zz - 1.0).abs() < tol
461            && self.xy.abs() < tol
462            && self.xz.abs() < tol
463            && self.yx.abs() < tol
464            && self.yz.abs() < tol
465            && self.zx.abs() < tol
466            && self.zy.abs() < tol
467    }
468
469    /// Linear interpolation
470    #[inline]
471    pub fn lerp(a: Self, b: Self, t: f64) -> Self {
472        let ot = 1.0 - t;
473        Self {
474            xx: ot * a.xx + t * b.xx,
475            xy: ot * a.xy + t * b.xy,
476            xz: ot * a.xz + t * b.xz,
477            yx: ot * a.yx + t * b.yx,
478            yy: ot * a.yy + t * b.yy,
479            yz: ot * a.yz + t * b.yz,
480            zx: ot * a.zx + t * b.zx,
481            zy: ot * a.zy + t * b.zy,
482            zz: ot * a.zz + t * b.zz,
483        }
484    }
485}
486
487// --- Conversions ---
488
489impl From<SphericalTensor> for Tensor {
490    #[inline]
491    fn from(st: SphericalTensor) -> Self {
492        Self {
493            xx: st.ii,
494            xy: 0.0,
495            xz: 0.0,
496            yx: 0.0,
497            yy: st.ii,
498            yz: 0.0,
499            zx: 0.0,
500            zy: 0.0,
501            zz: st.ii,
502        }
503    }
504}
505
506impl From<SymmTensor> for Tensor {
507    #[inline]
508    fn from(st: SymmTensor) -> Self {
509        Self {
510            xx: st.xx,
511            xy: st.xy,
512            xz: st.xz,
513            yx: st.xy,
514            yy: st.yy,
515            yz: st.yz,
516            zx: st.xz,
517            zy: st.yz,
518            zz: st.zz,
519        }
520    }
521}
522
523/// Hodge dual of a Vector as a skew-symmetric Tensor. C++ `operator*(Vector)`.
524pub fn hodge_dual_of_vec(v: Vector3) -> Tensor {
525    Tensor {
526        xx: 0.0,
527        xy: -v.z,
528        xz: v.y,
529        yx: v.z,
530        yy: 0.0,
531        yz: -v.x,
532        zx: -v.y,
533        zy: v.x,
534        zz: 0.0,
535    }
536}
537
538// --- Arithmetic operators (Tensor × Tensor, Tensor × scalar) ---
539
540impl Neg for Tensor {
541    type Output = Self;
542    #[inline]
543    fn neg(self) -> Self {
544        Self {
545            xx: -self.xx,
546            xy: -self.xy,
547            xz: -self.xz,
548            yx: -self.yx,
549            yy: -self.yy,
550            yz: -self.yz,
551            zx: -self.zx,
552            zy: -self.zy,
553            zz: -self.zz,
554        }
555    }
556}
557
558impl Add for Tensor {
559    type Output = Self;
560    #[inline]
561    fn add(self, r: Self) -> Self {
562        Self {
563            xx: self.xx + r.xx,
564            xy: self.xy + r.xy,
565            xz: self.xz + r.xz,
566            yx: self.yx + r.yx,
567            yy: self.yy + r.yy,
568            yz: self.yz + r.yz,
569            zx: self.zx + r.zx,
570            zy: self.zy + r.zy,
571            zz: self.zz + r.zz,
572        }
573    }
574}
575
576impl Sub for Tensor {
577    type Output = Self;
578    #[inline]
579    fn sub(self, r: Self) -> Self {
580        Self {
581            xx: self.xx - r.xx,
582            xy: self.xy - r.xy,
583            xz: self.xz - r.xz,
584            yx: self.yx - r.yx,
585            yy: self.yy - r.yy,
586            yz: self.yz - r.yz,
587            zx: self.zx - r.zx,
588            zy: self.zy - r.zy,
589            zz: self.zz - r.zz,
590        }
591    }
592}
593
594impl Mul<f64> for Tensor {
595    type Output = Self;
596    #[inline]
597    fn mul(self, s: f64) -> Self {
598        Self {
599            xx: self.xx * s,
600            xy: self.xy * s,
601            xz: self.xz * s,
602            yx: self.yx * s,
603            yy: self.yy * s,
604            yz: self.yz * s,
605            zx: self.zx * s,
606            zy: self.zy * s,
607            zz: self.zz * s,
608        }
609    }
610}
611
612impl Mul<Tensor> for f64 {
613    type Output = Tensor;
614    #[inline]
615    fn mul(self, t: Tensor) -> Tensor {
616        t * self
617    }
618}
619
620impl Div<f64> for Tensor {
621    type Output = Self;
622    #[inline]
623    fn div(self, s: f64) -> Self {
624        Self {
625            xx: self.xx / s,
626            xy: self.xy / s,
627            xz: self.xz / s,
628            yx: self.yx / s,
629            yy: self.yy / s,
630            yz: self.yz / s,
631            zx: self.zx / s,
632            zy: self.zy / s,
633            zz: self.zz / s,
634        }
635    }
636}
637
638impl AddAssign for Tensor {
639    #[inline]
640    fn add_assign(&mut self, r: Self) {
641        *self = *self + r;
642    }
643}
644
645impl SubAssign for Tensor {
646    #[inline]
647    fn sub_assign(&mut self, r: Self) {
648        *self = *self - r;
649    }
650}
651
652impl MulAssign<f64> for Tensor {
653    #[inline]
654    fn mul_assign(&mut self, s: f64) {
655        *self = *self * s;
656    }
657}
658
659// --- Mixed Tensor ± SphericalTensor ---
660
661impl Add<SphericalTensor> for Tensor {
662    type Output = Self;
663    #[inline]
664    fn add(self, st: SphericalTensor) -> Self {
665        Self {
666            xx: self.xx + st.ii,
667            xy: self.xy,
668            xz: self.xz,
669            yx: self.yx,
670            yy: self.yy + st.ii,
671            yz: self.yz,
672            zx: self.zx,
673            zy: self.zy,
674            zz: self.zz + st.ii,
675        }
676    }
677}
678
679impl Add<Tensor> for SphericalTensor {
680    type Output = Tensor;
681    #[inline]
682    fn add(self, t: Tensor) -> Tensor {
683        t + self
684    }
685}
686
687impl Sub<SphericalTensor> for Tensor {
688    type Output = Self;
689    #[inline]
690    fn sub(self, st: SphericalTensor) -> Self {
691        Self {
692            xx: self.xx - st.ii,
693            xy: self.xy,
694            xz: self.xz,
695            yx: self.yx,
696            yy: self.yy - st.ii,
697            yz: self.yz,
698            zx: self.zx,
699            zy: self.zy,
700            zz: self.zz - st.ii,
701        }
702    }
703}
704
705impl Sub<Tensor> for SphericalTensor {
706    type Output = Tensor;
707    #[inline]
708    fn sub(self, t: Tensor) -> Tensor {
709        Tensor {
710            xx: self.ii - t.xx,
711            xy: -t.xy,
712            xz: -t.xz,
713            yx: -t.yx,
714            yy: self.ii - t.yy,
715            yz: -t.yz,
716            zx: -t.zx,
717            zy: -t.zy,
718            zz: self.ii - t.zz,
719        }
720    }
721}
722
723// --- Mixed Tensor ± SymmTensor ---
724
725impl Add<SymmTensor> for Tensor {
726    type Output = Self;
727    #[inline]
728    fn add(self, st: SymmTensor) -> Self {
729        Self {
730            xx: self.xx + st.xx,
731            xy: self.xy + st.xy,
732            xz: self.xz + st.xz,
733            yx: self.yx + st.xy,
734            yy: self.yy + st.yy,
735            yz: self.yz + st.yz,
736            zx: self.zx + st.xz,
737            zy: self.zy + st.yz,
738            zz: self.zz + st.zz,
739        }
740    }
741}
742
743impl Add<Tensor> for SymmTensor {
744    type Output = Tensor;
745    #[inline]
746    fn add(self, t: Tensor) -> Tensor {
747        t + self
748    }
749}
750
751impl Sub<SymmTensor> for Tensor {
752    type Output = Self;
753    #[inline]
754    fn sub(self, st: SymmTensor) -> Self {
755        Self {
756            xx: self.xx - st.xx,
757            xy: self.xy - st.xy,
758            xz: self.xz - st.xz,
759            yx: self.yx - st.xy,
760            yy: self.yy - st.yy,
761            yz: self.yz - st.yz,
762            zx: self.zx - st.xz,
763            zy: self.zy - st.yz,
764            zz: self.zz - st.zz,
765        }
766    }
767}
768
769impl Sub<Tensor> for SymmTensor {
770    type Output = Tensor;
771    #[inline]
772    fn sub(self, t: Tensor) -> Tensor {
773        Tensor {
774            xx: self.xx - t.xx,
775            xy: self.xy - t.xy,
776            xz: self.xz - t.xz,
777            yx: self.xy - t.yx,
778            yy: self.yy - t.yy,
779            yz: self.yz - t.yz,
780            zx: self.xz - t.zx,
781            zy: self.yz - t.zy,
782            zz: self.zz - t.zz,
783        }
784    }
785}
786
787// SphericalTensor & Tensor → Tensor (scalar multiply of all elements)
788impl Mul<Tensor> for SphericalTensor {
789    type Output = Tensor;
790    #[inline]
791    fn mul(self, t: Tensor) -> Tensor {
792        t * self.ii
793    }
794}
795
796impl Mul<SphericalTensor> for Tensor {
797    type Output = Self;
798    #[inline]
799    fn mul(self, st: SphericalTensor) -> Self {
800        self * st.ii
801    }
802}
803
804// SymmTensor & Tensor → Tensor (matrix multiply)
805impl Mul<Tensor> for SymmTensor {
806    type Output = Tensor;
807    #[inline]
808    fn mul(self, t: Tensor) -> Tensor {
809        Tensor {
810            xx: self.xx * t.xx + self.xy * t.yx + self.xz * t.zx,
811            xy: self.xx * t.xy + self.xy * t.yy + self.xz * t.zy,
812            xz: self.xx * t.xz + self.xy * t.yz + self.xz * t.zz,
813            yx: self.xy * t.xx + self.yy * t.yx + self.yz * t.zx,
814            yy: self.xy * t.xy + self.yy * t.yy + self.yz * t.zy,
815            yz: self.xy * t.xz + self.yy * t.yz + self.yz * t.zz,
816            zx: self.xz * t.xx + self.yz * t.yx + self.zz * t.zx,
817            zy: self.xz * t.xy + self.yz * t.yy + self.zz * t.zy,
818            zz: self.xz * t.xz + self.yz * t.yz + self.zz * t.zz,
819        }
820    }
821}
822
823impl Mul<SymmTensor> for Tensor {
824    type Output = Self;
825    #[inline]
826    fn mul(self, st: SymmTensor) -> Self {
827        Self {
828            xx: self.xx * st.xx + self.xy * st.xy + self.xz * st.xz,
829            xy: self.xx * st.xy + self.xy * st.yy + self.xz * st.yz,
830            xz: self.xx * st.xz + self.xy * st.yz + self.xz * st.zz,
831            yx: self.yx * st.xx + self.yy * st.xy + self.yz * st.xz,
832            yy: self.yx * st.xy + self.yy * st.yy + self.yz * st.yz,
833            yz: self.yx * st.xz + self.yy * st.yz + self.yz * st.zz,
834            zx: self.zx * st.xx + self.zy * st.xy + self.zz * st.xz,
835            zy: self.zx * st.xy + self.zy * st.yy + self.zz * st.yz,
836            zz: self.zx * st.xz + self.zy * st.yz + self.zz * st.zz,
837        }
838    }
839}
840
841// SymmTensor & SymmTensor → Tensor (matrix multiply, result is NOT symmetric)
842impl Mul<SymmTensor> for SymmTensor {
843    type Output = Tensor;
844    #[inline]
845    fn mul(self, st: SymmTensor) -> Tensor {
846        Tensor {
847            xx: self.xx * st.xx + self.xy * st.xy + self.xz * st.xz,
848            xy: self.xx * st.xy + self.xy * st.yy + self.xz * st.yz,
849            xz: self.xx * st.xz + self.xy * st.yz + self.xz * st.zz,
850            yx: self.xy * st.xx + self.yy * st.xy + self.yz * st.xz,
851            yy: self.xy * st.xy + self.yy * st.yy + self.yz * st.yz,
852            yz: self.xy * st.xz + self.yy * st.yz + self.yz * st.zz,
853            zx: self.xz * st.xx + self.yz * st.xy + self.zz * st.xz,
854            zy: self.xz * st.xy + self.yz * st.yy + self.zz * st.yz,
855            zz: self.xz * st.xz + self.yz * st.yz + self.zz * st.zz,
856        }
857    }
858}
859
860// Outer product: Vector ⊗ Vector → Tensor. C++ `operator*(Vector, Vector)`.
861impl Mul<Vector3> for Vector3 {
862    type Output = Tensor;
863    #[inline]
864    fn mul(self, v: Vector3) -> Tensor {
865        Tensor {
866            xx: self.x * v.x,
867            xy: self.x * v.y,
868            xz: self.x * v.z,
869            yx: self.y * v.x,
870            yy: self.y * v.y,
871            yz: self.y * v.z,
872            zx: self.z * v.x,
873            zy: self.z * v.y,
874            zz: self.z * v.z,
875        }
876    }
877}
878
879// --- Free functions ---
880
881/// Trace tr(T) = xx + yy + zz.
882#[inline]
883pub fn tr(t: Tensor) -> f64 {
884    t.tr()
885}
886
887/// Determinant det(T).
888#[inline]
889pub fn det(t: Tensor) -> f64 {
890    t.det()
891}
892
893/// Inverse T⁻¹ (panics in debug builds if singular).
894#[inline]
895pub fn inv(t: Tensor) -> Tensor {
896    t.inv()
897}
898
899/// Symmetric part `0.5*(T + Tᵀ)`.
900#[inline]
901pub fn symm(t: Tensor) -> SymmTensor {
902    t.symm()
903}
904
905/// Twice the symmetric part `T + Tᵀ`.
906#[inline]
907pub fn two_symm(t: Tensor) -> SymmTensor {
908    t.two_symm()
909}
910
911/// Skew-symmetric part `0.5*(T - Tᵀ)`.
912#[inline]
913pub fn skew(t: Tensor) -> Tensor {
914    t.skew()
915}
916
917/// Deviatoric part `T - (tr/3)*I`.
918#[inline]
919pub fn dev(t: Tensor) -> Tensor {
920    t.dev()
921}
922
923/// Two-thirds deviatoric `T - (2*tr/3)*I`.
924#[inline]
925pub fn dev2(t: Tensor) -> Tensor {
926    t.dev2()
927}
928
929/// Deviatoric of the symmetric part `symm(T) - (tr/3)*I`.
930#[inline]
931pub fn dev_symm(t: Tensor) -> SymmTensor {
932    t.dev_symm()
933}
934
935/// Deviatoric of twice the symmetric part `twoSymm(T) - (2*tr/3)*I`.
936#[inline]
937pub fn dev_two_symm(t: Tensor) -> SymmTensor {
938    t.dev_two_symm()
939}
940
941/// Linear interpolation `(1-t)*a + t*b` between two tensors.
942#[inline]
943pub fn lerp(a: Tensor, b: Tensor, t: f64) -> Tensor {
944    Tensor::lerp(a, b, t)
945}
946
947/// Outer product v ⊗ w. Same as `v * w` but as a named function.
948#[inline]
949pub fn outer(v: Vector3, w: Vector3) -> Tensor {
950    v * w
951}
952
953#[cfg(test)]
954mod tests {
955    use super::*;
956
957    #[test]
958    fn identity_properties() {
959        let i = Tensor::IDENTITY;
960        assert_eq!(i.tr(), 3.0);
961        assert_eq!(i.det(), 1.0);
962        let inv = i.inv();
963        assert!((inv.xx - 1.0).abs() < 1e-14);
964    }
965
966    #[test]
967    fn mat_mul_identity() {
968        let i = Tensor::IDENTITY;
969        let t = Tensor::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
970        assert_eq!(i.mat_mul(t), t);
971        assert_eq!(t.mat_mul(i), t);
972    }
973
974    #[test]
975    fn mat_vec_identity() {
976        let v = Vector3::new(1.0, 2.0, 3.0);
977        assert_eq!(Tensor::IDENTITY.mat_vec(v), v);
978    }
979
980    #[test]
981    fn transpose() {
982        let t = Tensor::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
983        let tt = t.transpose();
984        assert_eq!(tt.xy, t.yx);
985        assert_eq!(tt.yx, t.xy);
986        assert_eq!(tt.xz, t.zx);
987    }
988
989    #[test]
990    fn symm_and_skew_sum_to_original() {
991        let t = Tensor::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
992        let s: Tensor = Tensor::from(t.symm());
993        let sk = t.skew();
994        let sum = s + sk;
995        assert!((sum.xx - t.xx).abs() < 1e-14);
996        assert!((sum.xy - t.xy).abs() < 1e-14);
997        assert!((sum.yz - t.yz).abs() < 1e-14);
998    }
999
1000    #[test]
1001    fn outer_product() {
1002        let x = Vector3::X;
1003        let y = Vector3::Y;
1004        let xy = x * y;
1005        assert_eq!(xy.xx, 0.0);
1006        assert_eq!(xy.xy, 1.0);
1007        assert_eq!(xy.yx, 0.0);
1008    }
1009
1010    #[test]
1011    fn dev_traceless() {
1012        let t = Tensor::new(4.0, 1.0, 0.0, 1.0, 2.0, 0.0, 0.0, 0.0, 3.0);
1013        assert!((t.dev().tr()).abs() < 1e-14);
1014    }
1015
1016    #[test]
1017    fn inv_roundtrip_is_identity() {
1018        // T · T⁻¹ should equal I (up to floating-point noise).
1019        let t = Tensor::new(2.0, 1.0, 0.0, 1.0, 3.0, 1.0, 0.0, 1.0, 2.0);
1020        assert!(
1021            t.mat_mul(t.inv()).is_identity(1e-12),
1022            "T · T⁻¹ is not identity"
1023        );
1024    }
1025
1026    #[test]
1027    fn double_inner_is_symmetric() {
1028        // T1:T2 == T2:T1 (double inner product / Frobenius inner product)
1029        let t1 = Tensor::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
1030        let t2 = Tensor::new(2.0, 0.0, 1.0, -1.0, 3.0, 0.0, 2.0, 1.0, 4.0);
1031        let diff = (t1.double_inner(t2) - t2.double_inner(t1)).abs();
1032        assert!(diff < 1e-12, "double_inner asymmetry = {diff:.3e}");
1033    }
1034
1035    #[test]
1036    fn dev2_is_two_thirds_trace_removed() {
1037        // dev2(T) = T − (2/3)·tr(T)·I (OpenFOAM convention)
1038        let t = Tensor::new(6.0, 1.0, 0.0, 1.0, 3.0, 0.0, 0.0, 0.0, 3.0);
1039        let tr = t.tr(); // 12.0
1040        let d = t.dev2();
1041        assert!((d.xx - (t.xx - 2.0 / 3.0 * tr)).abs() < 1e-14);
1042        assert!((d.yy - (t.yy - 2.0 / 3.0 * tr)).abs() < 1e-14);
1043        assert!((d.zz - (t.zz - 2.0 / 3.0 * tr)).abs() < 1e-14);
1044        assert!((d.xy - t.xy).abs() < 1e-14);
1045        // dev2 is NOT trace-free: tr(dev2) = tr(T) - 2*tr(T) = -tr(T)
1046        assert!(
1047            (d.tr() - (-tr)).abs() < 1e-14,
1048            "tr(dev2) = {:.3e}, expected {:.3e}",
1049            d.tr(),
1050            -tr
1051        );
1052    }
1053}
outram_foam_basic_lib/primitives/tensor.rs

outram_foam_basic_lib/primitives/
tensor.rs
