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

outram_foam_basic_lib/primitives/
tensor.rs
