russell_tensor/
spectral2.rs

1use super::{vec_dyad_vec, Mandel, Tensor2, SQRT_2, SQRT_3, SQRT_6};
2use crate::StrError;
3use russell_lab::{mat_eigen_sym_jacobi, Matrix, Vector};
4
5/// Holds the spectral representation of a symmetric second-order tensor
6pub struct Spectral2 {
7    /// The mandel representation
8    mandel: Mandel,
9
10    /// Holds the eigenvalues; dim = 3
11    pub lambda: Vector,
12
13    /// Holds the eigenprojectors; set of 3 symmetric Tensor2 (dim = 6 or 4)
14    pub projectors: Vec<Tensor2>,
15}
16
17impl Spectral2 {
18    /// Returns a new instance
19    ///
20    /// **Note:** Must call [Spectral2::compose] to calculate `lambda` and `projectors`.
21    pub fn new(two_dim: bool) -> Self {
22        let mandel = if two_dim {
23            Mandel::Symmetric2D
24        } else {
25            Mandel::Symmetric
26        };
27        Spectral2 {
28            mandel,
29            lambda: Vector::new(3),
30            projectors: vec![Tensor2::new(mandel), Tensor2::new(mandel), Tensor2::new(mandel)],
31        }
32    }
33
34    /// Performs the spectral decomposition of a symmetric second-order tensor
35    ///
36    /// # Results
37    ///
38    /// The results are available in [Spectral2::lambda] and [Spectral2::projectors].
39    pub fn decompose(&mut self, tt: &Tensor2) -> Result<(), StrError> {
40        if tt.mandel != self.mandel {
41            return Err("the mandel representation is incompatible");
42        }
43        let dim = tt.vec.dim();
44        if dim == 4 {
45            // eigenvalues and eigenvectors
46            let (t22, mut a) = tt.as_matrix_2d();
47            let mut l = Vector::new(2);
48            let mut v = Matrix::new(2, 2);
49            mat_eigen_sym_jacobi(&mut l, &mut v, &mut a)?;
50            self.lambda[0] = l[0];
51            self.lambda[1] = l[1];
52            self.lambda[2] = t22;
53
54            // extract eigenvectors
55            let u0 = Vector::from(&[v.get(0, 0), v.get(1, 0)]);
56            let u1 = Vector::from(&[v.get(0, 1), v.get(1, 1)]);
57
58            // compute eigenprojectors
59            vec_dyad_vec(&mut self.projectors[0], 1.0, &u0, &u0).unwrap();
60            vec_dyad_vec(&mut self.projectors[1], 1.0, &u1, &u1).unwrap();
61            self.projectors[2].clear();
62            self.projectors[2].vec[2] = 1.0;
63        } else {
64            // eigenvalues and eigenvectors
65            let mut a = tt.as_matrix();
66            let mut v = Matrix::new(3, 3);
67            mat_eigen_sym_jacobi(&mut self.lambda, &mut v, &mut a)?;
68
69            // extract eigenvectors
70            let u0 = Vector::from(&[v.get(0, 0), v.get(1, 0), v.get(2, 0)]);
71            let u1 = Vector::from(&[v.get(0, 1), v.get(1, 1), v.get(2, 1)]);
72            let u2 = Vector::from(&[v.get(0, 2), v.get(1, 2), v.get(2, 2)]);
73
74            // compute eigenprojectors
75            vec_dyad_vec(&mut self.projectors[0], 1.0, &u0, &u0).unwrap();
76            vec_dyad_vec(&mut self.projectors[1], 1.0, &u1, &u1).unwrap();
77            vec_dyad_vec(&mut self.projectors[2], 1.0, &u2, &u2).unwrap();
78        }
79        Ok(())
80    }
81
82    /// Composes a new tensor from the eigenprojectors and diagonal values (lambda)
83    pub fn compose(&self, composed: &mut Tensor2, lambda: &Vector) -> Result<(), StrError> {
84        if composed.mandel != self.mandel {
85            return Err("the mandel representation is incompatible");
86        }
87        if lambda.dim() != 3 {
88            return Err("lambda.dim must be equal to 3");
89        }
90        let n = self.projectors[0].vec.dim();
91        for i in 0..n {
92            composed.vec[i] = lambda[0] * self.projectors[0].vec[i]
93                + lambda[1] * self.projectors[1].vec[i]
94                + lambda[2] * self.projectors[2].vec[i];
95        }
96        Ok(())
97    }
98
99    /// Calculates the octahedral basis on the principal values space
100    ///
101    /// Returns `(λ_star_1, λ_star_2, λ_star_3)`
102    pub fn octahedral_basis(&self) -> (f64, f64, f64) {
103        let (s1, s2, s3) = (self.lambda[0], self.lambda[1], self.lambda[2]);
104        let ls1 = (2.0 * s1 - s2 - s3) / SQRT_6;
105        let ls2 = (s1 + s2 + s3) / SQRT_3;
106        let ls3 = (s3 - s2) / SQRT_2;
107        (ls1, ls2, ls3)
108    }
109}
110
111////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
112
113#[cfg(test)]
114mod tests {
115    use super::Spectral2;
116    use crate::{Mandel, SampleTensor2, SamplesTensor2, Tensor2, SQRT_3, SQRT_3_BY_2};
117    use russell_lab::{approx_eq, mat_approx_eq, vec_approx_eq, Matrix, Vector};
118
119    #[test]
120    fn decompose_captures_errors() {
121        let mut spec = Spectral2::new(false);
122        let tt = Tensor2::new(Mandel::General);
123        assert_eq!(
124            spec.decompose(&tt).err(),
125            Some("the mandel representation is incompatible")
126        );
127    }
128
129    #[test]
130    fn compose_capture_errors() {
131        let spec = Spectral2::new(false);
132        let mut tt = Tensor2::new(Mandel::Symmetric2D);
133        assert_eq!(
134            spec.compose(&mut tt, &spec.lambda).err(),
135            Some("the mandel representation is incompatible")
136        );
137        let mut tt = Tensor2::new(Mandel::Symmetric);
138        let lambda = Vector::new(1);
139        assert_eq!(
140            spec.compose(&mut tt, &lambda).err(),
141            Some("lambda.dim must be equal to 3")
142        );
143    }
144
145    fn check(spec: &mut Spectral2, sample: &SampleTensor2, tol_lambda: f64, tol_proj: f64, tol_spectral: f64) {
146        let correct_lambda = sample.eigenvalues.unwrap();
147        let correct_projectors = sample.eigenprojectors.unwrap();
148
149        // perform spectral decomposition of symmetric matrix
150        let mandel = spec.projectors[0].mandel;
151        let tt = Tensor2::from_matrix(&sample.matrix, mandel).unwrap();
152        spec.decompose(&tt).unwrap();
153
154        // print results
155        // println!("a =\n{}", tt.as_matrix());
156        // println!("λ = {}, {}, {}", spec.lambda[0], spec.lambda[1], spec.lambda[2]);
157        // println!("P0 =\n{}", spec.projectors[0].as_matrix());
158        // println!("P1 =\n{}", spec.projectors[1].as_matrix());
159        // println!("P2 =\n{}", spec.projectors[2].as_matrix());
160
161        // check eigenvalues
162        vec_approx_eq(&spec.lambda, &correct_lambda, tol_lambda);
163
164        // check eigenprojectors
165        let pp0 = spec.projectors[0].as_matrix();
166        let pp1 = spec.projectors[1].as_matrix();
167        let pp2 = spec.projectors[2].as_matrix();
168        let correct0 = Matrix::from(&correct_projectors[0]);
169        let correct1 = Matrix::from(&correct_projectors[1]);
170        let correct2 = Matrix::from(&correct_projectors[2]);
171        mat_approx_eq(&correct0, &pp0, tol_proj);
172        mat_approx_eq(&correct1, &pp1, tol_proj);
173        mat_approx_eq(&correct2, &pp2, tol_proj);
174
175        // compose
176        let mut tt_new = Tensor2::new(mandel);
177        spec.compose(&mut tt_new, &spec.lambda).unwrap();
178        let a_new = tt_new.as_matrix();
179        let a = Matrix::from(&sample.matrix);
180        mat_approx_eq(&a, &a_new, tol_spectral);
181    }
182
183    #[test]
184    fn decompose_and_compose_work_3d() {
185        let mut spec = Spectral2::new(false);
186        check(&mut spec, &SamplesTensor2::TENSOR_O, 1e-15, 1e-15, 1e-15);
187        check(&mut spec, &SamplesTensor2::TENSOR_I, 1e-15, 1e-15, 1e-15);
188        check(&mut spec, &SamplesTensor2::TENSOR_X, 1e-15, 1e-15, 1e-15);
189        check(&mut spec, &SamplesTensor2::TENSOR_Y, 1e-13, 1e-15, 1e-15);
190        check(&mut spec, &SamplesTensor2::TENSOR_Z, 1e-14, 1e-15, 1e-15);
191        check(&mut spec, &SamplesTensor2::TENSOR_U, 1e-13, 1e-15, 1e-14);
192        check(&mut spec, &SamplesTensor2::TENSOR_S, 1e-13, 1e-14, 1e-14);
193    }
194
195    #[test]
196    fn decompose_and_compose_work_2d() {
197        let mut spec = Spectral2::new(true);
198        check(&mut spec, &SamplesTensor2::TENSOR_O, 1e-15, 1e-15, 1e-15);
199        check(&mut spec, &SamplesTensor2::TENSOR_I, 1e-15, 1e-15, 1e-15);
200        check(&mut spec, &SamplesTensor2::TENSOR_X, 1e-15, 1e-15, 1e-15);
201        check(&mut spec, &SamplesTensor2::TENSOR_Y, 1e-13, 1e-15, 1e-15);
202        check(&mut spec, &SamplesTensor2::TENSOR_Z, 1e-14, 1e-15, 1e-15);
203    }
204
205    #[test]
206    fn octahedral_basis_works() {
207        // the following data corresponds to sigma_m = 1 and sigma_d = 3
208        #[rustfmt::skip]
209        let principal_stresses_and_lode = [
210            ( 3.0          ,  0.0          ,  0.0          ,  1.0 ),
211            ( 0.0          ,  3.0          ,  0.0          ,  1.0 ),
212            ( 0.0          ,  0.0          ,  3.0          ,  1.0 ),
213            ( 1.0 + SQRT_3 ,  1.0 - SQRT_3 ,  1.0          ,  0.0 ),
214            ( 1.0 + SQRT_3 ,  1.0          ,  1.0 - SQRT_3 ,  0.0 ),
215            ( 1.0          ,  1.0 + SQRT_3 ,  1.0 - SQRT_3 ,  0.0 ),
216            ( 1.0 - SQRT_3 ,  1.0 + SQRT_3 ,  1.0          ,  0.0 ),
217            ( 1.0          ,  1.0 - SQRT_3 ,  1.0 + SQRT_3 ,  0.0 ),
218            ( 1.0 - SQRT_3 ,  1.0          ,  1.0 + SQRT_3 ,  0.0 ),
219            ( 2.0          , -1.0          ,  2.0          , -1.0 ),
220            ( 2.0          ,  2.0          , -1.0          , -1.0 ),
221            (-1.0          ,  2.0          ,  2.0          , -1.0 ),
222        ];
223        let two_dim = true;
224        let mut spec = Spectral2::new(two_dim);
225        let mut tt = Tensor2::new_sym(two_dim);
226        for (sigma_1, sigma_2, sigma_3, lode_correct) in &principal_stresses_and_lode {
227            tt.vec[0] = *sigma_1;
228            tt.vec[1] = *sigma_2;
229            tt.vec[2] = *sigma_3;
230            spec.decompose(&tt).unwrap();
231            let (ls1, ls2, ls3) = spec.octahedral_basis();
232            let radius = f64::sqrt(ls3 * ls3 + ls1 * ls1);
233            let distance = ls2;
234            approx_eq(distance / SQRT_3, 1.0, 1e-15);
235            approx_eq(radius * SQRT_3_BY_2, 3.0, 1e-15);
236            if radius > 0.0 {
237                let cos_theta = ls1 / radius;
238                let lode = 4.0 * f64::powf(cos_theta, 3.0) - 3.0 * cos_theta;
239                approx_eq(lode, *lode_correct, 1e-15);
240            }
241        }
242    }
243}
russell_tensor/spectral2.rs

russell_tensor/
spectral2.rs
