russell_tensor/operations_mix1.rs
1use super::{Tensor2, Tensor4};
2use russell_lab::{mat_vec_mul, mat_vec_mul_update, vec_mat_mul, vec_outer, vec_outer_update};
3
4#[allow(unused)]
5use crate::Mandel; // for documentation
6
7/// Performs the dyadic product between two Tensor2 resulting a Tensor4
8///
9/// ```text
10/// D = α a ⊗ b
11/// ```
12///
13/// With orthonormal Cartesian components:
14///
15/// ```text
16/// Dᵢⱼₖₗ = α aᵢⱼ bₖₗ
17/// ```
18///
19/// Or, in Mandel basis:
20///
21/// ```text
22/// Dₘₙ = α aₘ bₙ
23/// ```
24///
25/// # Output
26///
27/// * `dd` -- the tensor `D`; with the same [Mandel] as `a` and `b`
28///
29/// # Input
30///
31/// * `a` -- first tensor; with the same [Mandel] as `b` and `dd`
32/// * `b` -- second tensor; with the same [Mandel] as `a` and `dd`
33///
34/// # Panics
35///
36/// A panic will occur the tensors have different [Mandel]
37///
38/// # Examples
39///
40/// ```
41/// use russell_tensor::{t2_dyad_t2, Mandel, Tensor2, Tensor4, StrError};
42///
43/// fn main() -> Result<(), StrError> {
44/// let a = Tensor2::from_matrix(&[
45/// [ 1.0, 10.0, 0.0],
46/// [-2.0, -1.0, 0.0],
47/// [ 0.0, 0.0, 2.0],
48/// ], Mandel::General)?;
49///
50/// let b = Tensor2::from_matrix(&[
51/// [1.0, 4.0, 6.0],
52/// [7.0, 2.0, 5.0],
53/// [9.0, 8.0, 3.0],
54/// ], Mandel::General)?;
55///
56/// let mut dd = Tensor4::new(Mandel::General);
57/// t2_dyad_t2(&mut dd, 1.0, &a, &b);
58///
59/// assert_eq!(
60/// format!("{:.1}", dd.as_matrix()),
61/// "┌ ┐\n\
62/// │ 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 │\n\
63/// │ -1.0 -2.0 -3.0 -4.0 -5.0 -6.0 -7.0 -8.0 -9.0 │\n\
64/// │ 2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 18.0 │\n\
65/// │ 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 │\n\
66/// │ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 │\n\
67/// │ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 │\n\
68/// │ -2.0 -4.0 -6.0 -8.0 -10.0 -12.0 -14.0 -16.0 -18.0 │\n\
69/// │ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 │\n\
70/// │ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 │\n\
71/// └ ┘"
72/// );
73/// Ok(())
74/// }
75/// ```
76pub fn t2_dyad_t2(dd: &mut Tensor4, alpha: f64, a: &Tensor2, b: &Tensor2) {
77 assert_eq!(a.mandel, dd.mandel);
78 assert_eq!(b.mandel, dd.mandel);
79 vec_outer(&mut dd.mat, alpha, &a.vec, &b.vec).unwrap();
80}
81
82/// Performs the dyadic product between two Tensor2 resulting in a Tensor4 (with update)
83///
84/// Computes:
85///
86/// ```text
87/// D += α a ⊗ b
88/// ```
89///
90/// With orthonormal Cartesian components:
91///
92/// ```text
93/// Dᵢⱼₖₗ += α aᵢⱼ bₖₗ
94/// ```
95///
96/// Or, in Mandel basis:
97///
98/// ```text
99/// Dₘₙ += α aₘ bₙ
100/// ```
101///
102/// # Output
103///
104/// * `dd` -- the tensor `D`; with the same [Mandel] as `a` and `b`
105///
106/// # Input
107///
108/// * `a` -- first tensor; with the same [Mandel] as `b` and `dd`
109/// * `b` -- second tensor; with the same [Mandel] as `a` and `dd`
110///
111/// # Panics
112///
113/// A panic will occur the tensors have different [Mandel]
114///
115/// # Examples
116///
117/// ```
118/// use russell_lab::Matrix;
119/// use russell_tensor::{t2_dyad_t2_update, Mandel, StrError, Tensor2, Tensor4};
120///
121/// fn main() -> Result<(), StrError> {
122/// #[rustfmt::skip]
123/// let a = Tensor2::from_matrix(&[
124/// [ 1.0, 10.0, 0.0],
125/// [ 2.0, 1.0, 0.0],
126/// [ 0.0, 0.0, 2.0],
127/// ], Mandel::General)?;
128///
129/// #[rustfmt::skip]
130/// let b = Tensor2::from_matrix(&[
131/// [1.0, 4.0, 6.0],
132/// [7.0, 2.0, 5.0],
133/// [9.0, 8.0, 3.0],
134/// ], Mandel::General)?;
135///
136/// let mat = Matrix::filled(9, 9, 0.5);
137/// let mut dd = Tensor4::from_matrix(&mat, Mandel::General)?;
138/// t2_dyad_t2_update(&mut dd, 1.0, &a, &b);
139///
140/// assert_eq!(
141/// format!("{:.1}", dd.as_matrix()),
142/// "┌ ┐\n\
143/// │ 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 │\n\
144/// │ 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 │\n\
145/// │ 2.5 4.5 6.5 8.5 10.5 12.5 14.5 16.5 18.5 │\n\
146/// │ 10.5 20.5 30.5 40.5 50.5 60.5 70.5 80.5 90.5 │\n\
147/// │ 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 │\n\
148/// │ 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 │\n\
149/// │ 2.5 4.5 6.5 8.5 10.5 12.5 14.5 16.5 18.5 │\n\
150/// │ 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 │\n\
151/// │ 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 │\n\
152/// └ ┘"
153/// );
154/// Ok(())
155/// }
156/// ```
157pub fn t2_dyad_t2_update(dd: &mut Tensor4, alpha: f64, a: &Tensor2, b: &Tensor2) {
158 assert_eq!(a.mandel, dd.mandel);
159 assert_eq!(b.mandel, dd.mandel);
160 vec_outer_update(&mut dd.mat, alpha, &a.vec, &b.vec).unwrap();
161}
162
163/// Performs the double-dot (ddot) operation between a Tensor4 and a Tensor2
164///
165/// ```text
166/// b = α D : a
167/// ```
168///
169/// With orthonormal Cartesian components:
170///
171/// ```text
172/// bᵢⱼ = α Σ Σ Dᵢⱼₖₗ aₖₗ
173/// k l
174/// ```
175///
176/// Or, in Mandel basis:
177///
178/// ```text
179/// bₘ = α Σ Dₘₙ aₙ
180/// n
181/// ```
182///
183/// # Output
184///
185/// * `b` -- the resulting second-order tensor; with the same [Mandel] as `a` and `dd`
186///
187/// # Input
188///
189/// * `alpha` -- the scalar multiplier
190/// * `dd` -- the fourth-order tensor; with the same [Mandel] as `a` and `b`
191/// * `a` -- the input second-order tensor; with the same [Mandel] as `b` and `dd`
192///
193/// # Panics
194///
195/// A panic will occur the tensors have different [Mandel]
196///
197/// # Examples
198///
199/// ```
200/// use russell_tensor::{t4_ddot_t2, Mandel, Tensor2, Tensor4, StrError};
201///
202/// fn main() -> Result<(), StrError> {
203/// let dd = Tensor4::from_matrix(&[
204/// [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0],
205/// [ -1.0, -2.0, -3.0, -4.0, -5.0, -6.0, -7.0, -8.0, -9.0],
206/// [ 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0],
207/// [ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0],
208/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
209/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
210/// [ -2.0, -4.0, -6.0, -8.0,-10.0,-12.0,-14.0,-16.0,-18.0],
211/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
212/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
213/// ], Mandel::General)?;
214///
215/// let a = Tensor2::from_matrix(&[
216/// [1.0, 4.0, 6.0],
217/// [7.0, 2.0, 5.0],
218/// [9.0, 8.0, 3.0],
219/// ], Mandel::General)?;
220///
221/// let mut b = Tensor2::new(Mandel::General);
222/// t4_ddot_t2(&mut b, 1.0, &dd, &a);
223///
224/// assert_eq!(
225/// format!("{:.1}", b.as_matrix()),
226/// "┌ ┐\n\
227/// │ 285.0 2850.0 0.0 │\n\
228/// │ -570.0 -285.0 0.0 │\n\
229/// │ 0.0 0.0 570.0 │\n\
230/// └ ┘"
231/// );
232/// Ok(())
233/// }
234/// ```
235pub fn t4_ddot_t2(b: &mut Tensor2, alpha: f64, dd: &Tensor4, a: &Tensor2) {
236 assert_eq!(a.mandel, dd.mandel);
237 assert_eq!(b.mandel, dd.mandel);
238 mat_vec_mul(&mut b.vec, alpha, &dd.mat, &a.vec).unwrap();
239}
240
241/// Performs the double-dot (ddot) operation between a Tensor4 and a Tensor2 with update
242///
243/// Computes:
244///
245/// ```text
246/// b = α D : a + β b
247/// ```
248///
249/// With orthonormal Cartesian components:
250///
251/// ```text
252/// bᵢⱼ = α Σ Σ Dᵢⱼₖₗ aₖₗ + β bᵢⱼ
253/// k l
254/// ```
255///
256/// Or, in Mandel basis:
257///
258/// ```text
259/// bₘ = α Σ Dₘₙ aₙ + β bₘ
260/// n
261/// ```
262///
263/// # Output
264///
265/// * `b` -- the resulting second-order tensor; with the same [Mandel] as `a` and `dd`
266///
267/// # Input
268///
269/// * `alpha` -- the scalar multiplier
270/// * `a` -- the input second-order tensor; with the same [Mandel] as `b` and `dd`
271/// * `dd` -- the fourth-order tensor; with the same [Mandel] as `a` and `b`
272/// * `beta` -- the other scalar multiplier
273///
274/// # Panics
275///
276/// A panic will occur the tensors have different [Mandel]
277///
278/// # Examples
279///
280/// ```
281/// use russell_tensor::{t4_ddot_t2_update, Mandel, Tensor2, Tensor4, StrError};
282///
283/// fn main() -> Result<(), StrError> {
284/// let dd = Tensor4::from_matrix(&[
285/// [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0],
286/// [ -1.0, -2.0, -3.0, -4.0, -5.0, -6.0, -7.0, -8.0, -9.0],
287/// [ 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0],
288/// [ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0],
289/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
290/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
291/// [ -2.0, -4.0, -6.0, -8.0,-10.0,-12.0,-14.0,-16.0,-18.0],
292/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
293/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
294/// ], Mandel::General)?;
295///
296/// let a = Tensor2::from_matrix(&[
297/// [1.0, 4.0, 6.0],
298/// [7.0, 2.0, 5.0],
299/// [9.0, 8.0, 3.0],
300/// ], Mandel::General)?;
301///
302/// let mut b = Tensor2::from_matrix(&[
303/// [1.0, 0.0, 0.0],
304/// [0.0, 1.0, 0.0],
305/// [0.0, 0.0, 1.0],
306/// ], Mandel::General)?;
307/// t4_ddot_t2_update(&mut b, 1.0, &dd, &a, 1000.0);
308///
309/// assert_eq!(
310/// format!("{:.1}", b.as_matrix()),
311/// "┌ ┐\n\
312/// │ 1285.0 2850.0 0.0 │\n\
313/// │ -570.0 715.0 0.0 │\n\
314/// │ 0.0 0.0 1570.0 │\n\
315/// └ ┘"
316/// );
317/// Ok(())
318/// }
319/// ```
320pub fn t4_ddot_t2_update(b: &mut Tensor2, alpha: f64, dd: &Tensor4, a: &Tensor2, beta: f64) {
321 assert_eq!(a.mandel, dd.mandel);
322 assert_eq!(b.mandel, dd.mandel);
323 mat_vec_mul_update(&mut b.vec, alpha, &dd.mat, &a.vec, beta).unwrap();
324}
325
326/// Performs the double-dot (ddot) operation between a Tensor2 and a Tensor4
327///
328/// Computes:
329///
330/// ```text
331/// b = α a : D
332/// ```
333///
334/// With orthonormal Cartesian components:
335///
336/// ```text
337/// bₖₗ = α Σ Σ aᵢⱼ Dᵢⱼₖₗ
338/// i j
339/// ```
340///
341/// # Output
342///
343/// * `b` -- the resulting second-order tensor; with the same [Mandel] as `a` and `dd`
344///
345/// # Input
346///
347/// * `alpha` -- the scalar multiplier
348/// * `a` -- the input second-order tensor; with the same [Mandel] as `b` and `dd`
349/// * `dd` -- the fourth-order tensor; with the same [Mandel] as `a` and `b`
350///
351/// # Panics
352///
353/// A panic will occur the tensors have different [Mandel]
354///
355/// # Examples
356///
357/// ```
358/// use russell_tensor::{t2_ddot_t4, Mandel, Tensor2, Tensor4, StrError};
359///
360/// fn main() -> Result<(), StrError> {
361/// let a = Tensor2::from_matrix(&[
362/// [1.0, 4.0, 6.0],
363/// [7.0, 2.0, 5.0],
364/// [9.0, 8.0, 3.0],
365/// ], Mandel::General)?;
366///
367/// let dd = Tensor4::from_matrix(&[
368/// [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0],
369/// [ -1.0, -2.0, -3.0, -4.0, -5.0, -6.0, -7.0, -8.0, -9.0],
370/// [ 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0],
371/// [ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0],
372/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
373/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
374/// [ -2.0, -4.0, -6.0, -8.0,-10.0,-12.0,-14.0,-16.0,-18.0],
375/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
376/// [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
377/// ], Mandel::General)?;
378///
379/// let mut b = Tensor2::new(Mandel::General);
380/// t2_ddot_t4(&mut b, 1.0, &a, &dd);
381///
382/// assert_eq!(
383/// format!("{:.1}", b.as_matrix()),
384/// "┌ ┐\n\
385/// │ 31.0 124.0 186.0 │\n\
386/// │ 217.0 62.0 155.0 │\n\
387/// │ 279.0 248.0 93.0 │\n\
388/// └ ┘"
389/// );
390/// Ok(())
391/// }
392/// ```
393pub fn t2_ddot_t4(b: &mut Tensor2, alpha: f64, a: &Tensor2, dd: &Tensor4) {
394 assert_eq!(a.mandel, dd.mandel);
395 assert_eq!(b.mandel, dd.mandel);
396 vec_mat_mul(&mut b.vec, alpha, &a.vec, &dd.mat).unwrap();
397}
398
399/// Computes Tensor2 double-dot Tensor4 double-dot Tensor2
400///
401/// Computes:
402///
403/// ```text
404/// s = a : D : b
405/// ```
406///
407/// With orthonormal Cartesian components:
408///
409/// ```text
410/// s = Σ Σ Σ Σ aᵢⱼ Dᵢⱼₖₗ bₖₗ
411/// i j k l
412/// ```
413///
414/// Or, in Mandel basis:
415///
416/// ```text
417/// s = Σ Σ aₘ Dₘₙ bₙ
418/// m n
419/// ```
420///
421/// Note: the Lagrange multiplier in Plasticity needs this operation.
422///
423/// # Input
424///
425/// * `a` -- the first second-order tensor; with the same [Mandel] as `b` and `dd`
426/// * `dd` -- the fourth-order tensor; with the same [Mandel] as `a` and `b`
427/// * `b` -- the second second-order tensor; with the same [Mandel] as `a` and `dd`
428///
429/// # Output
430///
431/// Returns the scalar results.
432///
433/// # Panics
434///
435/// A panic will occur the tensors have different [Mandel]
436pub fn t2_ddot_t4_ddot_t2(a: &Tensor2, dd: &Tensor4, b: &Tensor2) -> f64 {
437 assert_eq!(a.mandel, dd.mandel);
438 assert_eq!(b.mandel, dd.mandel);
439 let dim = a.vec.dim();
440 let mut s = 0.0;
441 for m in 0..dim {
442 for n in 0..dim {
443 s += a.vec[m] * dd.mat.get(m, n) * b.vec[n];
444 }
445 }
446 s
447}
448
449/// Computes Tensor4 double-dot Tensor2 dyadic Tensor2 double-dot Tensor4
450///
451/// Computes:
452///
453/// ```text
454/// E = α D + β (D : a) ⊗ (b : D)
455/// ```
456///
457/// With orthonormal Cartesian components:
458///
459/// ```text
460/// Eᵢⱼₖₗ = α Dᵢⱼₖₗ + β Σ Σ Σ Σ (Dᵢⱼₛₜ aₛₜ) (bₒₚ Dₒₚₖₗ)
461/// s t o p
462/// ```
463///
464/// Or, in Mandel basis:
465///
466/// ```text
467/// Eₘₙ = α Dₘₙ + β Σ Σ (Dₘₐ aₐ) (bₑ Dₑₙ)
468/// a e
469/// ```
470///
471/// Note: the elastoplastic modulus in Plasticity needs this operation.
472///
473/// # Output
474///
475/// * `ee` -- the resulting fourth-order tensor; with the same [Mandel] as the other tensors
476///
477/// # Input
478///
479/// * `alpha` -- the first scalar multiplier
480/// * `dd` -- the fourth-order tensor; with the same [Mandel] as the other tensors
481/// * `beta` -- the second scalar multiplier
482/// * `a` -- the first second-order tensor; with the same [Mandel] as the other tensors
483/// * `b` -- the second second-order tensor; with the same [Mandel] as the other tensors
484///
485/// # Panics
486///
487/// A panic will occur the tensors have different [Mandel]
488pub fn t4_ddot_t2_dyad_t2_ddot_t4(ee: &mut Tensor4, alpha: f64, dd: &Tensor4, beta: f64, a: &Tensor2, b: &Tensor2) {
489 assert_eq!(a.mandel, dd.mandel);
490 assert_eq!(b.mandel, dd.mandel);
491 assert_eq!(ee.mandel, dd.mandel);
492 let dim = a.vec.dim();
493 for m in 0..dim {
494 for n in 0..dim {
495 ee.mat.set(m, n, alpha * dd.mat.get(m, n));
496 for p in 0..dim {
497 for q in 0..dim {
498 ee.mat
499 .add(m, n, beta * dd.mat.get(m, p) * a.vec[p] * b.vec[q] * dd.mat.get(q, n));
500 }
501 }
502 }
503 }
504}
505
506////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
507
508#[cfg(test)]
509mod tests {
510 use super::*;
511 use crate::{Mandel, SamplesTensor4, MN_TO_IJKL};
512 use russell_lab::{approx_eq, mat_approx_eq, Matrix};
513
514 #[test]
515 #[should_panic]
516 fn t2_dyad_t2_panics_on_different_mandel1() {
517 let a = Tensor2::new(Mandel::Symmetric2D);
518 let b = Tensor2::new(Mandel::Symmetric); // wrong; it must be Symmetric2D
519 let mut dd = Tensor4::new(Mandel::Symmetric2D);
520 t2_dyad_t2(&mut dd, 1.0, &a, &b);
521 }
522
523 #[test]
524 #[should_panic]
525 fn t2_dyad_t2_panics_on_different_mandel2() {
526 let a = Tensor2::new(Mandel::Symmetric2D);
527 let b = Tensor2::new(Mandel::Symmetric2D);
528 let mut dd = Tensor4::new(Mandel::Symmetric); // wrong; it must be Symmetric2D
529 t2_dyad_t2(&mut dd, 1.0, &a, &b);
530 }
531
532 #[test]
533 fn t2_dyad_t2_works() {
534 // general dyad general
535 #[rustfmt::skip]
536 let a = Tensor2::from_matrix(&[
537 [1.0, 2.0, 3.0],
538 [4.0, 5.0, 6.0],
539 [7.0, 8.0, 9.0],
540 ], Mandel::General).unwrap();
541 #[rustfmt::skip]
542 let b = Tensor2::from_matrix(&[
543 [0.5, 0.5, 0.5],
544 [0.5, 0.5, 0.5],
545 [0.5, 0.5, 0.5],
546 ], Mandel::General).unwrap();
547 let mut dd = Tensor4::new(Mandel::General);
548 t2_dyad_t2(&mut dd, 2.0, &a, &b);
549 let mat = dd.as_matrix();
550 assert_eq!(
551 format!("{:.1}", mat),
552 "┌ ┐\n\
553 │ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 │\n\
554 │ 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 │\n\
555 │ 9.0 9.0 9.0 9.0 9.0 9.0 9.0 9.0 9.0 │\n\
556 │ 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 │\n\
557 │ 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 │\n\
558 │ 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 │\n\
559 │ 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 │\n\
560 │ 8.0 8.0 8.0 8.0 8.0 8.0 8.0 8.0 8.0 │\n\
561 │ 7.0 7.0 7.0 7.0 7.0 7.0 7.0 7.0 7.0 │\n\
562 └ ┘"
563 );
564
565 // sym-3D dyad general
566 #[rustfmt::skip]
567 let a = Tensor2::from_matrix(&[
568 [1.0, 2.0, 3.0],
569 [2.0, 5.0, 6.0],
570 [3.0, 6.0, 9.0],
571 ], Mandel::Symmetric).unwrap();
572 #[rustfmt::skip]
573 let b = Tensor2::from_matrix(&[
574 [0.5, 0.5, 0.5],
575 [0.5, 0.5, 0.5],
576 [0.5, 0.5, 0.5],
577 ], Mandel::Symmetric).unwrap();
578 let mut dd = Tensor4::new(Mandel::Symmetric);
579 t2_dyad_t2(&mut dd, 2.0, &a, &b);
580 let mat = dd.as_matrix();
581 assert_eq!(
582 format!("{:.1}", mat),
583 "┌ ┐\n\
584 │ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 │\n\
585 │ 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 │\n\
586 │ 9.0 9.0 9.0 9.0 9.0 9.0 9.0 9.0 9.0 │\n\
587 │ 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 │\n\
588 │ 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 │\n\
589 │ 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 │\n\
590 │ 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 │\n\
591 │ 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 │\n\
592 │ 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 │\n\
593 └ ┘"
594 );
595
596 // sym-2D dyad sym-2D
597 #[rustfmt::skip]
598 let a = Tensor2::from_matrix(&[
599 [1.0, 2.0, 0.0],
600 [2.0, 5.0, 0.0],
601 [0.0, 0.0, 9.0],
602 ], Mandel::Symmetric2D).unwrap();
603 #[rustfmt::skip]
604 let b = Tensor2::from_matrix(&[
605 [0.5, 0.5, 0.0],
606 [0.5, 0.5, 0.0],
607 [0.0, 0.0, 0.5],
608 ], Mandel::Symmetric2D).unwrap();
609 let mut dd = Tensor4::new(Mandel::Symmetric2D);
610 t2_dyad_t2(&mut dd, 2.0, &a, &b);
611 let mat = dd.as_matrix();
612 assert_eq!(
613 format!("{:.1}", mat),
614 "┌ ┐\n\
615 │ 1.0 1.0 1.0 1.0 0.0 0.0 1.0 0.0 0.0 │\n\
616 │ 5.0 5.0 5.0 5.0 0.0 0.0 5.0 0.0 0.0 │\n\
617 │ 9.0 9.0 9.0 9.0 0.0 0.0 9.0 0.0 0.0 │\n\
618 │ 2.0 2.0 2.0 2.0 0.0 0.0 2.0 0.0 0.0 │\n\
619 │ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 │\n\
620 │ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 │\n\
621 │ 2.0 2.0 2.0 2.0 0.0 0.0 2.0 0.0 0.0 │\n\
622 │ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 │\n\
623 │ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 │\n\
624 └ ┘"
625 );
626 }
627
628 #[test]
629 #[should_panic]
630 fn t2_dyad_t2_update_panics_on_different_mandel1() {
631 let a = Tensor2::new(Mandel::Symmetric2D);
632 let b = Tensor2::new(Mandel::Symmetric); // wrong; it must be Symmetric2D
633 let mut dd = Tensor4::new(Mandel::Symmetric2D);
634 t2_dyad_t2_update(&mut dd, 1.0, &a, &b);
635 }
636
637 #[test]
638 #[should_panic]
639 fn t2_dyad_t2_update_panics_on_different_mandel2() {
640 let a = Tensor2::new(Mandel::Symmetric2D);
641 let b = Tensor2::new(Mandel::Symmetric2D);
642 let mut dd = Tensor4::new(Mandel::Symmetric); // wrong; it must be Symmetric2D
643 t2_dyad_t2_update(&mut dd, 1.0, &a, &b);
644 }
645
646 #[test]
647 fn t2_dyad_t2_update_works() {
648 // general dyad general
649 #[rustfmt::skip]
650 let a = Tensor2::from_matrix(&[
651 [1.0, 2.0, 3.0],
652 [4.0, 5.0, 6.0],
653 [7.0, 8.0, 9.0],
654 ], Mandel::General).unwrap();
655 #[rustfmt::skip]
656 let b = Tensor2::from_matrix(&[
657 [0.5, 0.5, 0.5],
658 [0.5, 0.5, 0.5],
659 [0.5, 0.5, 0.5],
660 ], Mandel::General).unwrap();
661 let mat = Matrix::filled(9, 9, 0.1);
662 let mut dd = Tensor4::from_matrix(&mat, Mandel::General).unwrap();
663 t2_dyad_t2_update(&mut dd, 2.0, &a, &b);
664 let mat = dd.as_matrix();
665 let correct = "┌ ┐\n\
666 │ 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 │\n\
667 │ 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.1 │\n\
668 │ 9.1 9.1 9.1 9.1 9.1 9.1 9.1 9.1 9.1 │\n\
669 │ 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 2.1 │\n\
670 │ 6.1 6.1 6.1 6.1 6.1 6.1 6.1 6.1 6.1 │\n\
671 │ 3.1 3.1 3.1 3.1 3.1 3.1 3.1 3.1 3.1 │\n\
672 │ 4.1 4.1 4.1 4.1 4.1 4.1 4.1 4.1 4.1 │\n\
673 │ 8.1 8.1 8.1 8.1 8.1 8.1 8.1 8.1 8.1 │\n\
674 │ 7.1 7.1 7.1 7.1 7.1 7.1 7.1 7.1 7.1 │\n\
675 └ ┘";
676 assert_eq!(format!("{:.1}", mat), correct);
677 }
678
679 fn check_dyad(s: f64, a_ten: &Tensor2, b_ten: &Tensor2, dd_ten: &Tensor4, tol: f64) {
680 let a = a_ten.as_matrix();
681 let b = b_ten.as_matrix();
682 let dd = dd_ten.as_matrix();
683 let mut correct = Matrix::new(9, 9);
684 for m in 0..9 {
685 for n in 0..9 {
686 let (i, j, k, l) = MN_TO_IJKL[m][n];
687 correct.set(m, n, s * a.get(i, j) * b.get(k, l));
688 }
689 }
690 mat_approx_eq(&dd, &correct, tol);
691 }
692
693 #[test]
694 fn t2_dyad_t2_works_extra() {
695 // general dyad general
696 #[rustfmt::skip]
697 let a = Tensor2::from_matrix(&[
698 [1.0, 2.0, 3.0],
699 [4.0, 5.0, 6.0],
700 [7.0, 8.0, 9.0],
701 ], Mandel::General).unwrap();
702 #[rustfmt::skip]
703 let b = Tensor2::from_matrix(&[
704 [9.0, 8.0, 7.0],
705 [6.0, 5.0, 4.0],
706 [3.0, 2.0, 1.0],
707 ], Mandel::General).unwrap();
708 let mut dd = Tensor4::new(Mandel::General);
709 t2_dyad_t2(&mut dd, 2.0, &a, &b);
710 let mat = dd.as_matrix();
711 // println!("{:.1}", mat);
712 let correct = Matrix::from(&[
713 [18.0, 10.0, 2.0, 16.0, 8.0, 14.0, 12.0, 4.0, 6.0],
714 [90.0, 50.0, 10.0, 80.0, 40.0, 70.0, 60.0, 20.0, 30.0],
715 [162.0, 90.0, 18.0, 144.0, 72.0, 126.0, 108.0, 36.0, 54.0],
716 [36.0, 20.0, 4.0, 32.0, 16.0, 28.0, 24.0, 8.0, 12.0],
717 [108.0, 60.0, 12.0, 96.0, 48.0, 84.0, 72.0, 24.0, 36.0],
718 [54.0, 30.0, 6.0, 48.0, 24.0, 42.0, 36.0, 12.0, 18.0],
719 [72.0, 40.0, 8.0, 64.0, 32.0, 56.0, 48.0, 16.0, 24.0],
720 [144.0, 80.0, 16.0, 128.0, 64.0, 112.0, 96.0, 32.0, 48.0],
721 [126.0, 70.0, 14.0, 112.0, 56.0, 98.0, 84.0, 28.0, 42.0],
722 ]);
723 mat_approx_eq(&mat, &correct, 1e-13);
724 check_dyad(2.0, &a, &b, &dd, 1e-13);
725
726 // symmetric dyad symmetric
727 #[rustfmt::skip]
728 let a = Tensor2::from_matrix(&[
729 [1.0, 4.0, 6.0],
730 [4.0, 2.0, 5.0],
731 [6.0, 5.0, 3.0],
732 ], Mandel::Symmetric).unwrap();
733 #[rustfmt::skip]
734 let b = Tensor2::from_matrix(&[
735 [3.0, 5.0, 6.0],
736 [5.0, 2.0, 4.0],
737 [6.0, 4.0, 1.0],
738 ], Mandel::Symmetric).unwrap();
739 let mut dd = Tensor4::new(Mandel::Symmetric);
740 t2_dyad_t2(&mut dd, 2.0, &a, &b);
741 let mat = dd.as_matrix();
742 // println!("{:.1}", mat);
743 let correct = Matrix::from(&[
744 [6.0, 4.0, 2.0, 10.0, 8.0, 12.0, 10.0, 8.0, 12.0],
745 [12.0, 8.0, 4.0, 20.0, 16.0, 24.0, 20.0, 16.0, 24.0],
746 [18.0, 12.0, 6.0, 30.0, 24.0, 36.0, 30.0, 24.0, 36.0],
747 [24.0, 16.0, 8.0, 40.0, 32.0, 48.0, 40.0, 32.0, 48.0],
748 [30.0, 20.0, 10.0, 50.0, 40.0, 60.0, 50.0, 40.0, 60.0],
749 [36.0, 24.0, 12.0, 60.0, 48.0, 72.0, 60.0, 48.0, 72.0],
750 [24.0, 16.0, 8.0, 40.0, 32.0, 48.0, 40.0, 32.0, 48.0],
751 [30.0, 20.0, 10.0, 50.0, 40.0, 60.0, 50.0, 40.0, 60.0],
752 [36.0, 24.0, 12.0, 60.0, 48.0, 72.0, 60.0, 48.0, 72.0],
753 ]);
754 mat_approx_eq(&mat, &correct, 1e-13);
755 check_dyad(2.0, &a, &b, &dd, 1e-13);
756
757 // symmetric 2D dyad symmetric 2D
758 #[rustfmt::skip]
759 let a = Tensor2::from_matrix(&[
760 [1.0, 4.0, 0.0],
761 [4.0, 2.0, 0.0],
762 [0.0, 0.0, 3.0],
763 ], Mandel::Symmetric2D).unwrap();
764 #[rustfmt::skip]
765 let b = Tensor2::from_matrix(&[
766 [3.0, 4.0, 0.0],
767 [4.0, 2.0, 0.0],
768 [0.0, 0.0, 1.0],
769 ], Mandel::Symmetric2D).unwrap();
770 let mut dd = Tensor4::new(Mandel::Symmetric2D);
771 t2_dyad_t2(&mut dd, 2.0, &a, &b);
772 let mat = dd.as_matrix();
773 // println!("{:.1}", mat);
774 let correct = Matrix::from(&[
775 [6.0, 4.0, 2.0, 8.0, 0.0, 0.0, 8.0, 0.0, 0.0],
776 [12.0, 8.0, 4.0, 16.0, 0.0, 0.0, 16.0, 0.0, 0.0],
777 [18.0, 12.0, 6.0, 24.0, 0.0, 0.0, 24.0, 0.0, 0.0],
778 [24.0, 16.0, 8.0, 32.0, 0.0, 0.0, 32.0, 0.0, 0.0],
779 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
780 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
781 [24.0, 16.0, 8.0, 32.0, 0.0, 0.0, 32.0, 0.0, 0.0],
782 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
783 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
784 ]);
785 mat_approx_eq(&mat, &correct, 1e-14);
786 check_dyad(2.0, &a, &b, &dd, 1e-15);
787 }
788
789 #[test]
790 #[should_panic]
791 fn t4_ddot_t2_panics_on_different_mandel1() {
792 let a = Tensor2::new(Mandel::Symmetric); // wrong; it must be the same as `dd`
793 let mut b = Tensor2::new(Mandel::Symmetric2D);
794 let dd = Tensor4::new(Mandel::Symmetric2D);
795 t4_ddot_t2(&mut b, 1.0, &dd, &a);
796 }
797
798 #[test]
799 #[should_panic]
800 fn t4_ddot_t2_panics_on_different_mandel2() {
801 let a = Tensor2::new(Mandel::Symmetric2D);
802 let mut b = Tensor2::new(Mandel::Symmetric); // wrong; it must be the same as `dd`
803 let dd = Tensor4::new(Mandel::Symmetric2D);
804 t4_ddot_t2(&mut b, 1.0, &dd, &a);
805 }
806
807 #[test]
808 fn t4_ddot_t2_works() {
809 let dd = Tensor4::from_matrix(&SamplesTensor4::SYM_2D_SAMPLE1_STD_MATRIX, Mandel::Symmetric2D).unwrap();
810 #[rustfmt::skip]
811 let a = Tensor2::from_matrix(&[
812 [-1.0, -2.0, 0.0],
813 [-2.0, 2.0, 0.0],
814 [ 0.0, 0.0, -3.0]], Mandel::Symmetric2D).unwrap();
815 let mut b = Tensor2::new(Mandel::Symmetric2D);
816 t4_ddot_t2(&mut b, 1.0, &dd, &a);
817 let out = b.as_matrix();
818 assert_eq!(
819 format!("{:.1}", out),
820 "┌ ┐\n\
821 │ -46.0 -154.0 0.0 │\n\
822 │ -154.0 -64.0 0.0 │\n\
823 │ 0.0 0.0 -82.0 │\n\
824 └ ┘"
825 );
826 }
827
828 #[test]
829 #[should_panic]
830 fn t4_ddot_t2_update_panics_on_different_mandel1() {
831 let a = Tensor2::new(Mandel::Symmetric); // wrong; it must be the same as `dd`
832 let mut b = Tensor2::new(Mandel::Symmetric2D);
833 let dd = Tensor4::new(Mandel::Symmetric2D);
834 t4_ddot_t2_update(&mut b, 1.0, &dd, &a, 1.0);
835 }
836
837 #[test]
838 #[should_panic]
839 fn t4_ddot_update_t2_panics_on_different_mandel2() {
840 let a = Tensor2::new(Mandel::Symmetric2D);
841 let mut b = Tensor2::new(Mandel::Symmetric); // wrong; it must be the same as `dd`
842 let dd = Tensor4::new(Mandel::Symmetric2D);
843 t4_ddot_t2_update(&mut b, 1.0, &dd, &a, 1.0);
844 }
845
846 #[test]
847 fn t4_ddot_t2_update_works() {
848 let dd = Tensor4::from_matrix(&SamplesTensor4::SYM_2D_SAMPLE1_STD_MATRIX, Mandel::Symmetric2D).unwrap();
849 #[rustfmt::skip]
850 let a = Tensor2::from_matrix(&[
851 [-1.0, -2.0, 0.0],
852 [-2.0, 2.0, 0.0],
853 [ 0.0, 0.0, -3.0],
854 ], Mandel::Symmetric2D).unwrap();
855 #[rustfmt::skip]
856 let mut b = Tensor2::from_matrix(&[
857 [-1000.0, -1000.0, 0.0],
858 [-1000.0, -1000.0, 0.0],
859 [ 0.0, 0.0, -1000.0],
860 ], Mandel::Symmetric2D).unwrap();
861 t4_ddot_t2_update(&mut b, 1.0, &dd, &a, 2.0);
862 let out = b.as_matrix();
863 assert_eq!(
864 format!("{:.1}", out),
865 "┌ ┐\n\
866 │ -2046.0 -2154.0 0.0 │\n\
867 │ -2154.0 -2064.0 0.0 │\n\
868 │ 0.0 0.0 -2082.0 │\n\
869 └ ┘"
870 );
871 }
872
873 #[test]
874 #[should_panic]
875 fn t2_ddot_t4_panics_on_different_mandel1() {
876 let a = Tensor2::new(Mandel::Symmetric); // wrong; it must be the same as `dd`
877 let mut b = Tensor2::new(Mandel::Symmetric2D);
878 let dd = Tensor4::new(Mandel::Symmetric2D);
879 t2_ddot_t4(&mut b, 1.0, &a, &dd);
880 }
881
882 #[test]
883 #[should_panic]
884 fn t2_ddot_t4_panics_on_different_mandel2() {
885 let a = Tensor2::new(Mandel::Symmetric2D);
886 let mut b = Tensor2::new(Mandel::Symmetric); // wrong; it must be the same as `dd`
887 let dd = Tensor4::new(Mandel::Symmetric2D);
888 t2_ddot_t4(&mut b, 1.0, &a, &dd);
889 }
890
891 #[test]
892 fn t2_ddot_t4_works() {
893 let dd = Tensor4::from_matrix(&SamplesTensor4::SYM_2D_SAMPLE1_STD_MATRIX, Mandel::Symmetric2D).unwrap();
894 #[rustfmt::skip]
895 let a = Tensor2::from_matrix(&[
896 [-1.0, -2.0, 0.0],
897 [-2.0, 2.0, 0.0],
898 [ 0.0, 0.0, -3.0]], Mandel::Symmetric2D).unwrap();
899 let mut b = Tensor2::new(Mandel::Symmetric2D);
900 t2_ddot_t4(&mut b, 1.0, &a, &dd);
901 let out = b.as_matrix();
902 assert_eq!(
903 format!("{:.1}", out),
904 "┌ ┐\n\
905 │ -90.0 -144.0 0.0 │\n\
906 │ -144.0 -96.0 0.0 │\n\
907 │ 0.0 0.0 -102.0 │\n\
908 └ ┘"
909 );
910 }
911
912 #[test]
913 #[should_panic]
914 fn t2_ddot_t4_ddot_t2_panics_on_different_mandel1() {
915 let a = Tensor2::new(Mandel::Symmetric); // wrong; it must be the same as `dd`
916 let b = Tensor2::new(Mandel::Symmetric2D);
917 let dd = Tensor4::new(Mandel::Symmetric2D);
918 t2_ddot_t4_ddot_t2(&a, &dd, &b);
919 }
920
921 #[test]
922 #[should_panic]
923 fn t2_ddot_t4_ddot_t2_panics_on_different_mandel2() {
924 let a = Tensor2::new(Mandel::Symmetric2D);
925 let b = Tensor2::new(Mandel::Symmetric); // wrong; it must be the same as `dd`
926 let dd = Tensor4::new(Mandel::Symmetric2D);
927 t2_ddot_t4_ddot_t2(&a, &dd, &b);
928 }
929
930 #[test]
931 fn t2_ddot_t4_ddot_t2_works() {
932 #[rustfmt::skip]
933 let a = Tensor2::from_matrix(&[
934 [1.0, 2.0, 3.0],
935 [4.0, 5.0, 6.0],
936 [7.0, 8.0, 9.0],
937 ], Mandel::General).unwrap();
938 #[rustfmt::skip]
939 let b = Tensor2::from_matrix(&[
940 [9.0, 8.0, 7.0],
941 [6.0, 5.0, 4.0],
942 [3.0, 2.0, 1.0],
943 ], Mandel::General).unwrap();
944 let mat = Matrix::filled(9, 9, -1.0);
945 let dd = Tensor4::from_matrix(&mat, Mandel::General).unwrap();
946 let s = t2_ddot_t4_ddot_t2(&a, &dd, &b);
947 approx_eq(s, -2025.0, 1e-15);
948 }
949
950 #[test]
951 #[should_panic]
952 fn t4_ddot_t2_dyad_t2_ddot_t4_panics_on_different_mandel1() {
953 let a = Tensor2::new(Mandel::Symmetric); // wrong; it must be the same as `ee`
954 let b = Tensor2::new(Mandel::Symmetric2D);
955 let dd = Tensor4::new(Mandel::Symmetric2D);
956 let mut ee = Tensor4::new(Mandel::Symmetric2D);
957 t4_ddot_t2_dyad_t2_ddot_t4(&mut ee, 2.0, &dd, 3.0, &a, &b);
958 }
959
960 #[test]
961 #[should_panic]
962 fn t4_ddot_t2_dyad_t2_ddot_t4_panics_on_different_mandel2() {
963 let a = Tensor2::new(Mandel::Symmetric2D);
964 let b = Tensor2::new(Mandel::Symmetric); // wrong; it must be the same as `ee`
965 let dd = Tensor4::new(Mandel::Symmetric2D);
966 let mut ee = Tensor4::new(Mandel::Symmetric2D);
967 t4_ddot_t2_dyad_t2_ddot_t4(&mut ee, 2.0, &dd, 3.0, &a, &b);
968 }
969
970 #[test]
971 #[should_panic]
972 fn t4_ddot_t2_dyad_t2_ddot_t4_panics_on_different_mandel3() {
973 let a = Tensor2::new(Mandel::Symmetric2D);
974 let b = Tensor2::new(Mandel::Symmetric2D);
975 let dd = Tensor4::new(Mandel::Symmetric); // wrong; it must be the same as `ee`
976 let mut ee = Tensor4::new(Mandel::Symmetric2D);
977 t4_ddot_t2_dyad_t2_ddot_t4(&mut ee, 2.0, &dd, 3.0, &a, &b);
978 }
979
980 #[test]
981 fn t4_ddot_t2_dyad_t2_ddot_t4_works1() {
982 #[rustfmt::skip]
983 let a = Tensor2::from_matrix(&[
984 [1.0, 2.0, 3.0],
985 [4.0, 5.0, 6.0],
986 [7.0, 8.0, 9.0],
987 ], Mandel::General).unwrap();
988 #[rustfmt::skip]
989 let b = Tensor2::from_matrix(&[
990 [9.0, 8.0, 7.0],
991 [6.0, 5.0, 4.0],
992 [3.0, 2.0, 1.0],
993 ], Mandel::General).unwrap();
994 let mat = Matrix::filled(9, 9, -1.0);
995 let dd = Tensor4::from_matrix(&mat, Mandel::General).unwrap();
996 let mut ee = Tensor4::new(Mandel::General);
997 t4_ddot_t2_dyad_t2_ddot_t4(&mut ee, 2.0, &dd, 3.0, &a, &b);
998 let correct = [
999 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1000 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1001 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1002 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1003 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1004 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1005 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1006 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1007 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1008 ];
1009 mat_approx_eq(&ee.as_matrix(), &correct, 1e-15);
1010 }
1011
1012 /*
1013 #[test]
1014 fn t4_ddot_t2_dyad_t2_ddot_t4_works2() {
1015 #[rustfmt::skip]
1016 let a = Tensor2::from_matrix(&[
1017 [1.0, 2.0, 3.0],
1018 [2.0, 5.0, 6.0],
1019 [3.0, 6.0, 9.0],
1020 ], Mandel::Symmetric).unwrap();
1021 #[rustfmt::skip]
1022 let b = Tensor2::from_matrix(&[
1023 [1.0, 4.0, 5.0],
1024 [4.0, 2.0, 6.0],
1025 [5.0, 6.0, 3.0],
1026 ], Mandel::Symmetric).unwrap();
1027 let dd = Tensor4::from_matrix(
1028 &[
1029 [1.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1030 [2.0, 4.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1031 [3.0, 3.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1032 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1033 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1034 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1035 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1036 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1037 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1038 ],
1039 Mandel::Symmetric,
1040 )
1041 .unwrap();
1042 let mut ee = Tensor4::new(Mandel::Symmetric);
1043 t4_ddot_t2_dyad_t2_ddot_t4(&mut ee, 2.0, &dd, 3.0, &a, &b);
1044 println!("{}", ee.as_matrix());
1045 let correct = [
1046 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1047 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1048 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1049 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1050 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1051 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1052 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1053 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1054 [6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073., 6073.],
1055 ];
1056 mat_approx_eq(&ee.as_matrix(), &correct, 1e-15);
1057 }
1058 */
1059}