russell_tensor 3.2.0

/// sqrt(2) <https://oeis.org/A002193>
#[allow(clippy::approx_constant)]
pub const SQRT_2: f64 =
    1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157f64;

/// sqrt(3) <https://oeis.org/A002194>
pub const SQRT_3: f64 = 1.7320508075688772935274463415058723669428052538103806280558069794519330169088000370811461867572485756756261414154f64;

/// sqrt(6) <https://oeis.org/A010464>
pub const SQRT_6: f64 =
    2.44948974278317809819728407470589139196594748065667012843269256725096037745731502653985943310464023f64;

/// sqrt(2/3) <https://oeis.org/A157697>
pub const SQRT_2_BY_3: f64 =
    0.816496580927726032732428024901963797321982493552223376144230855750320125819105008846619811034880078272864f64;

/// sqrt(3/2) <https://oeis.org/A115754>
pub const SQRT_3_BY_2: f64 =
    1.22474487139158904909864203735294569598297374032833506421634628362548018872865751326992971655232011f64;

/// 1/3
pub const ONE_BY_3: f64 =
    0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333f64;

/// 2/3
pub const TWO_BY_3: f64 =
    0.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666f64;

/// Identifies an operation as an addition (update)
///
/// Example u += A . v
pub const ADD: u8 = 1;

/// Identifies an operation as a setting
///
/// Example u = A . v
pub const SET: u8 = 0;

/// Tolerance to avoid zero division with the J2 invariant
///
/// This constant is used, for instance, in the calculation of the Lode invariant using the following expression:
///
/// ```text
///        3 √3 J3
/// l = ─────────────
///     2 pow(J2,1.5)
/// ```
///
/// Note: pow(1e-9,1.5) = 3.16e-14
pub const TOL_J2: f64 = 1e-9;

/// Indicates that a set of three Tensor2 satisfy the eigenprojectors rules
///
/// The set must satisfy idempotency (x), orthogonality (y), and completeness (z).
/// Each satisfied rule is assigned the number `1`. Thus, `xyz` must be `111`. This
/// constant is padded with an extra number such as `7` to facilitate automation.
/// Therefore the value `7111` means that all eigenprojector rules are satisfied.
pub const OK_EIGENPROJ_RULES: usize = 7111;

/// Holds the value 1.0 (f64); used to give the following constants the f64 type
const ONE: f64 = 1.0;

/// Second-order identity tensor in Kelvin-Mandel basis (I)
///
/// ```text
/// Kelvin-Mandel vector:
///       ┌   ┐
///       │ 1 │
///       │ 1 │
///       │ 1 │
///       │ 0 │
/// [I] = │ 0 │
///       │ 0 │
///       │ 0 │
///       │ 0 │
///       │ 0 │
///       └   ┘
/// ```
pub const IDENTITY2: [f64; 9] = [ONE, ONE, ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];

/// Fourth-order identity tensor in Kelvin-Mandel basis (II)
///
/// **Note:** this tensor cannot be represented in reduced-dimension because it is not minor-symmetric.
///
/// ```text
/// Definition:
///        _
/// II = I ⊗ I
/// ```
///
/// ```text
/// Kelvin-Mandel matrix:
///        ┌                     ┐
///        │ 1 0 0  0 0 0  0 0 0 │
///        │ 0 1 0  0 0 0  0 0 0 │
///        │ 0 0 1  0 0 0  0 0 0 │
///        │ 0 0 0  1 0 0  0 0 0 │
/// [II] = │ 0 0 0  0 1 0  0 0 0 │
///        │ 0 0 0  0 0 1  0 0 0 │
///        │ 0 0 0  0 0 0  1 0 0 │
///        │ 0 0 0  0 0 0  0 1 0 │
///        │ 0 0 0  0 0 0  0 0 1 │
///        └                     ┘
/// ```
pub const IDENTITY4: [[f64; 9]; 9] = [
    [ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ONE],
];

/// Fourth-order transposition tensor in Kelvin-Mandel basis (TT)
///
/// **Note:** this tensor cannot be represented in reduced-dimension because it is not minor-symmetric.
///
/// ```text
/// Definition:
///
/// TT = I ⊗ I
///        ‾
/// ```
///
/// ```text
/// Kelvin-Mandel matrix:
///        ┌                        ┐
///        │ 1 0 0  0 0 0   0  0  0 │
///        │ 0 1 0  0 0 0   0  0  0 │
///        │ 0 0 1  0 0 0   0  0  0 │
///        │ 0 0 0  1 0 0   0  0  0 │
/// [TT] = │ 0 0 0  0 1 0   0  0  0 │
///        │ 0 0 0  0 0 1   0  0  0 │
///        │ 0 0 0  0 0 0  -1  0  0 │
///        │ 0 0 0  0 0 0   0 -1  0 │
///        │ 0 0 0  0 0 0   0  0 -1 │
///        └                        ┘
/// ```
pub const TRANSPOSITION: [[f64; 9]; 9] = [
    [ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -ONE, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -ONE, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -ONE],
];

/// Fourth-order trace-projection tensor (JJ)
///
/// Note: this tensor can be represented in reduced-dimension.
///
/// ```text
/// Definition:
///
/// JJ = I ⊗ I
/// ```
///
/// ```text
/// Kelvin-Mandel matrix:
///        ┌                     ┐
///        │ 1 1 1  0 0 0  0 0 0 │
///        │ 1 1 1  0 0 0  0 0 0 │
///        │ 1 1 1  0 0 0  0 0 0 │
///        │ 0 0 0  0 0 0  0 0 0 │
/// [JJ] = │ 0 0 0  0 0 0  0 0 0 │
///        │ 0 0 0  0 0 0  0 0 0 │
///        │ 0 0 0  0 0 0  0 0 0 │
///        │ 0 0 0  0 0 0  0 0 0 │
///        │ 0 0 0  0 0 0  0 0 0 │
///        └                     ┘
/// ```
pub const TRACE_PROJECTION: [[f64; 9]; 9] = [
    [ONE, ONE, ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [ONE, ONE, ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [ONE, ONE, ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
];

/// Fourth-order isotropic making projector (Piso)
///
/// Note: this tensor can be represented in reduced-dimension.
///
/// ```text
/// Definition:
///
/// Piso = ⅓ I ⊗ I = ⅓ JJ
/// ```
///
/// ```text
/// Kelvin-Mandel matrix:
///          ┌                     ┐
///          │ ⅓ ⅓ ⅓  0 0 0  0 0 0 │
///          │ ⅓ ⅓ ⅓  0 0 0  0 0 0 │
///          │ ⅓ ⅓ ⅓  0 0 0  0 0 0 │
///          │ 0 0 0  0 0 0  0 0 0 │
/// [Piso] = │ 0 0 0  0 0 0  0 0 0 │
///          │ 0 0 0  0 0 0  0 0 0 │
///          │ 0 0 0  0 0 0  0 0 0 │
///          │ 0 0 0  0 0 0  0 0 0 │
///          │ 0 0 0  0 0 0  0 0 0 │
///          └                     ┘
/// ```
pub const P_ISO: [[f64; 9]; 9] = [
    [ONE_BY_3, ONE_BY_3, ONE_BY_3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [ONE_BY_3, ONE_BY_3, ONE_BY_3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [ONE_BY_3, ONE_BY_3, ONE_BY_3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
];

/// Fourth-order symmetric making projector (Psym)
///
/// Note: this tensor can be represented in reduced-dimension.
///
/// ```text
/// Definition:
///             _
/// Psym = ½ (I ⊗ I + I ⊗ I) = ½ (II + TT) = ½ ssd(I)
///                     ‾
/// ```
///
/// ```text
/// Kelvin-Mandel matrix:
///          ┌                     ┐
///          │ 1 0 0  0 0 0  0 0 0 │
///          │ 0 1 0  0 0 0  0 0 0 │
///          │ 0 0 1  0 0 0  0 0 0 │
///          │ 0 0 0  1 0 0  0 0 0 │
/// [Psym] = │ 0 0 0  0 1 0  0 0 0 │
///          │ 0 0 0  0 0 1  0 0 0 │
///          │ 0 0 0  0 0 0  0 0 0 │
///          │ 0 0 0  0 0 0  0 0 0 │
///          │ 0 0 0  0 0 0  0 0 0 │
///          └                     ┘
/// ```
pub const P_SYM: [[f64; 9]; 9] = [
    [ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
];

/// Fourth-order skew making projector (Pskew)
///
/// **Note:** this tensor cannot be represented in reduced-dimension because it is not minor-symmetric.
///
/// ```text
/// Definition:
///              _
/// Pskew = ½ (I ⊗ I - I ⊗ I) = ½ (II - TT)
///                      ‾
/// ```
///
/// ```text
/// Kelvin-Mandel matrix:
///           ┌                     ┐
///           │ 0 0 0  0 0 0  0 0 0 │
///           │ 0 0 0  0 0 0  0 0 0 │
///           │ 0 0 0  0 0 0  0 0 0 │
///           │ 0 0 0  0 0 0  0 0 0 │
/// [Pskew] = │ 0 0 0  0 0 0  0 0 0 │
///           │ 0 0 0  0 0 0  0 0 0 │
///           │ 0 0 0  0 0 0  1 0 0 │
///           │ 0 0 0  0 0 0  0 1 0 │
///           │ 0 0 0  0 0 0  0 0 1 │
///           └                     ┘
/// ```
pub const P_SKEW: [[f64; 9]; 9] = [
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ONE],
];

/// Fourth-order deviatoric making projector (Pdev)
///
/// **Note:** this tensor cannot be represented in reduced-dimension because it is not minor-symmetric.
///
/// ```text
/// Definition:
///          _
/// Pdev = I ⊗ I - ⅓ I ⊗ I = II - Piso
/// ```
///
/// ```text
/// Kelvin-Mandel matrix:
///          ┌                        ┐
///          │  ⅔ -⅓ -⅓  0 0 0  0 0 0 │
///          │ -⅓  ⅔ -⅓  0 0 0  0 0 0 │
///          │ -⅓ -⅓  ⅔  0 0 0  0 0 0 │
///          │  0  0  0  1 0 0  0 0 0 │
/// [Pdev] = │  0  0  0  0 1 0  0 0 0 │
///          │  0  0  0  0 0 1  0 0 0 │
///          │  0  0  0  0 0 0  1 0 0 │
///          │  0  0  0  0 0 0  0 1 0 │
///          │  0  0  0  0 0 0  0 0 1 │
///          └                        ┘
/// ```
pub const P_DEV: [[f64; 9]; 9] = [
    [TWO_BY_3, -ONE_BY_3, -ONE_BY_3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-ONE_BY_3, TWO_BY_3, -ONE_BY_3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-ONE_BY_3, -ONE_BY_3, TWO_BY_3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ONE],
];

/// Fourth-order symmetric-deviatoric projector in Kelvin-Mandel basis
///
/// Note: this tensor can be represented in reduced-dimension.
///
/// ```text
/// Definition:
///                _
/// Psymdev = ½ (I ⊗ I + I ⊗ I) - ⅓ I ⊗ I = Psym - Piso
///                        ‾
/// ```
///
/// ```text
/// Kelvin-Mandel matrix:
///             ┌                        ┐
///             │  ⅔ -⅓ -⅓  0 0 0  0 0 0 │
///             │ -⅓  ⅔ -⅓  0 0 0  0 0 0 │
///             │ -⅓ -⅓  ⅔  0 0 0  0 0 0 │
///             │  0  0  0  1 0 0  0 0 0 │
/// [Psymdev] = │  0  0  0  0 1 0  0 0 0 │
///             │  0  0  0  0 0 1  0 0 0 │
///             │  0  0  0  0 0 0  0 0 0 │
///             │  0  0  0  0 0 0  0 0 0 │
///             │  0  0  0  0 0 0  0 0 0 │
///             └                        ┘
/// ```
pub const P_SYMDEV: [[f64; 9]; 9] = [
    [TWO_BY_3, -ONE_BY_3, -ONE_BY_3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-ONE_BY_3, TWO_BY_3, -ONE_BY_3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [-ONE_BY_3, -ONE_BY_3, TWO_BY_3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, ONE, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
];

// --- maps -------------------------------------------------------------------------------------------------------------

//
// --- Second-order tensors ---
//

/// Maps the m-th position in the vector representation to the index (i,j) of Tensor2
///
/// Diagonal goes first, then the upper off-diagonal entries (i < j), and, finally, the lower off-diagonal entries (i > j).
///
/// ```text
/// ┌   ┐    ┌    ┐
/// │ 0 │    │ 00 │
/// │ 1 │    │ 11 │
/// │ 2 │    │ 22 │    ┌          ┐
/// │ 3 │    │ 01 │    │ 00 01 02 │
/// │ 4 │ => │ 12 │ => │ 10 11 12 │
/// │ 5 │    │ 02 │    │ 20 21 22 │
/// │ 6 │    │ 10 │    └          ┘
/// │ 7 │    │ 21 │
/// │ 8 │    │ 20 │
/// └   ┘    └    ┘
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::M_TO_IJ;
/// assert_eq!(M_TO_IJ[3], (0,1));
/// ```
#[rustfmt::skip]
pub const M_TO_IJ: [(usize, usize); 9] = [
    // diagonal
    (0,0), // 0
    (1,1), // 1
    (2,2), // 2
    // upper-diagonal
    (0,1), // 3
    (1,2), // 4
    (0,2), // 5
    // lower-diagonal
    (1,0), // 6
    (2,1), // 7
    (2,0), // 8
];

/// Maps (i,j) of Tensor2 to the m-th position in the vector representation
///
/// Diagonal goes first, then the upper off-diagonal entries (i < j), and, finally, the lower off-diagonal entries (i > j).
///
/// ```text
///                 ┌    ┐    ┌   ┐
///                 │ 00 │    │ 0 │
///                 │ 11 │    │ 1 │
/// ┌          ┐    │ 22 │    │ 2 │
/// │ 00 01 02 │    │ 01 │    │ 3 │
/// │ 10 11 12 │ => │ 12 │ => │ 4 │
/// │ 20 21 22 │    │ 02 │    │ 5 │
/// └          ┘    │ 10 │    │ 6 │
///                 │ 21 │    │ 7 │
///                 │ 20 │    │ 8 │
///                 └    ┘    └   ┘
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::IJ_TO_M;
/// assert_eq!(IJ_TO_M[0][1], 3);
/// ```
#[rustfmt::skip]
pub const IJ_TO_M: [[usize; 3]; 3] = [
    [0, 3, 5],
    [6, 1, 4],
    [8, 7, 2],
];

/// Maps (i,j) of Tensor2 to the m-th position in the vector representation (symmetric version)
///
/// Diagonal goes first, then the upper diagonals, and, finally, the lower diagonals.
///
/// ```text
///                 ┌    ┐    ┌   ┐
/// ┌          ┐    │ 00 │    │ 0 │
/// │ 00 01 02 │    │ 11 │    │ 1 │
/// │ 01 11 12 │ => │ 22 │ => │ 2 │
/// │ 02 12 22 │    │ 01 │    │ 3 │
/// └          ┘    │ 12 │    │ 4 │
///                 │ 02 │    │ 5 │
///                 └    ┘    └   ┘
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::IJ_TO_M_SYM;
/// assert_eq!(IJ_TO_M_SYM[0][1], 3);
/// ```
#[rustfmt::skip]
pub const IJ_TO_M_SYM: [[usize; 3]; 3] = [
    [0, 3, 5],
    [3, 1, 4],
    [5, 4, 2],
];

//
// --- Third-order tensors ---
//

/// Maps the (m,n)-th position in the matrix representation to (i,j,k) of Tensor3 (Case A)
///
/// * `m` pairs with the first two indices `ij`
/// * `n` equals `k`
/// 
/// ```text
/// Case A:
///      00    01    02
///    -----------------
/// 0 │ 00_0  00_1  00_2
/// 1 │ 11_0  11_1  11_2
/// 2 │ 22_0  22_1  22_2
/// 3 │ 01_0  01_1  01_2
/// 4 │ 12_0  12_1  12_2
/// 5 │ 02_0  02_1  02_2
/// 6 │ 10_0  10_1  10_2
/// 7 │ 21_0  21_1  21_2
/// 8 │ 20_0  20_1  20_2
///    -----------------
///      80    81    82
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::MN_TO_IJK_CASE_A;
/// assert_eq!(MN_TO_IJK_CASE_A[3][2], (0,1,2));
/// ```
#[rustfmt::skip]
pub const MN_TO_IJK_CASE_A: [[(usize,usize,usize); 3]; 9] = [
    [(0,0,0), (0,0,1), (0,0,2)], // 00_0  00_1  00_2
    [(1,1,0), (1,1,1), (1,1,2)], // 11_0  11_1  11_2
    [(2,2,0), (2,2,1), (2,2,2)], // 22_0  22_1  22_2
    [(0,1,0), (0,1,1), (0,1,2)], // 01_0  01_1  01_2
    [(1,2,0), (1,2,1), (1,2,2)], // 12_0  12_1  12_2
    [(0,2,0), (0,2,1), (0,2,2)], // 02_0  02_1  02_2
    [(1,0,0), (1,0,1), (1,0,2)], // 10_0  10_1  10_2
    [(2,1,0), (2,1,1), (2,1,2)], // 21_0  21_1  21_2
    [(2,0,0), (2,0,1), (2,0,2)], // 20_0  20_1  20_2
];

/// Maps the (m,n)-th position in the matrix representation to (i,j,k) of Tensor3 (Case B)
///
/// * `m` equals `i`
/// * `n` pairs with the last two indices `jk`
/// 
/// ```text
/// Case B:
///      00   01   02   03   04   05   06   07   08
///    ---------------------------------------------
/// 0 │ 0_00 0_11 0_22 0_01 0_12 0_02 0_10 0_21 0_20
/// 1 │ 1_00 1_11 1_22 1_01 1_12 1_02 1_10 1_21 1_20
/// 2 │ 2_00 2_11 2_22 2_01 2_12 2_02 2_10 2_21 2_20
///    ---------------------------------------------
///      20   21   22   23   24   25   26   27   28
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::MN_TO_IJK_CASE_B;
/// assert_eq!(MN_TO_IJK_CASE_B[2][3], (2,0,1));
/// ```
#[rustfmt::skip]
pub const MN_TO_IJK_CASE_B: [[(usize,usize,usize); 9]; 3] = [
    [(0,0,0), (0,1,1), (0,2,2), (0,0,1), (0,1,2), (0,0,2), (0,1,0), (0,2,1), (0,2,0)], // 0_00 0_11 0_22 0_01 0_12 0_02 0_10 0_21 0_20
    [(1,0,0), (1,1,1), (1,2,2), (1,0,1), (1,1,2), (1,0,2), (1,1,0), (1,2,1), (1,2,0)], // 1_00 1_11 1_22 1_01 1_12 1_02 1_10 1_21 1_20
    [(2,0,0), (2,1,1), (2,2,2), (2,0,1), (2,1,2), (2,0,2), (2,1,0), (2,2,1), (2,2,0)], // 2_00 2_11 2_22 2_01 2_12 2_02 2_10 2_21 2_20
];

/// Maps (i,j,k) of Tensor3 to the (m,n)-th position in the matrix representation (Case A)
///
/// * `m` pairs with the first two indices `ij`
/// * `n` equals `k`
/// 
/// ```text
/// Case A:
///      00    01    02
///    -----------------
/// 0 │ 00_0  00_1  00_2
/// 1 │ 11_0  11_1  11_2
/// 2 │ 22_0  22_1  22_2
/// 3 │ 01_0  01_1  01_2
/// 4 │ 12_0  12_1  12_2
/// 5 │ 02_0  02_1  02_2
/// 6 │ 10_0  10_1  10_2
/// 7 │ 21_0  21_1  21_2
/// 8 │ 20_0  20_1  20_2
///    -----------------
///      80    81    82
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::IJK_TO_MN_CASE_A;
/// assert_eq!(IJK_TO_MN_CASE_A[0][1][2], (3, 2));
/// ```
#[rustfmt::skip]
pub const IJK_TO_MN_CASE_A: [[[(usize, usize); 3]; 3]; 3] = [
    [
        [(0,0), (0,1), (0,2)], // 00_0 00_1 00_2
        [(3,0), (3,1), (3,2)], // 01_0 01_1 01_2
        [(5,0), (5,1), (5,2)], // 02_0 02_1 02_2
    ],
    [
        [(6,0), (6,1), (6,2)], // 10_0 10_1 10_2
        [(1,0), (1,1), (1,2)], // 11_0 11_1 11_2
        [(4,0), (4,1), (4,2)], // 12_0 12_1 12_2
    ],
    [
        [(8,0), (8,1), (8,2)], // 20_0 20_1 20_2
        [(7,0), (7,1), (7,2)], // 21_0 21_1 21_2
        [(2,0), (2,1), (2,2)], // 22_0 22_1 22_2
    ],
];

/// Maps (i,j,k) of Tensor3 to the (m,n)-th position in the matrix representation (Case B)
///
/// * `m` equals `i`
/// * `n` pairs with the last two indices `jk`
/// 
/// ```text
/// Case B:
///      00   01   02   03   04   05   06   07   08
///    ---------------------------------------------
/// 0 │ 0_00 0_11 0_22 0_01 0_12 0_02 0_10 0_21 0_20
/// 1 │ 1_00 1_11 1_22 1_01 1_12 1_02 1_10 1_21 1_20
/// 2 │ 2_00 2_11 2_22 2_01 2_12 2_02 2_10 2_21 2_20
///    ---------------------------------------------
///      20   21   22   23   24   25   26   27   28
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::IJK_TO_MN_CASE_B;
/// assert_eq!(IJK_TO_MN_CASE_B[0][1][2], (0, 4));
/// ```
#[rustfmt::skip]
pub const IJK_TO_MN_CASE_B: [[[(usize, usize); 3]; 3]; 3] = [
    [
        [(0,0),(0,3),(0,5)], // 0_00 0_01 0_02
        [(0,6),(0,1),(0,4)], // 0_10 0_11 0_12
        [(0,8),(0,7),(0,2)], // 0_20 0_21 0_22
    ],
    [
        [(1,0),(1,3),(1,5)], // 1_00 1_01 1_02
        [(1,6),(1,1),(1,4)], // 1_10 1_11 1_12
        [(1,8),(1,7),(1,2)], // 1_20 1_21 1_22
    ],
    [
        [(2,0),(2,3),(2,5)], // 2_00 2_01 2_02
        [(2,6),(2,1),(2,4)], // 2_10 2_11 2_12
        [(2,8),(2,7),(2,2)], // 2_20 2_21 2_22
    ],
];

/// Maps (i,j,k) of Tensor3 to the (m,n)-th position in the matrix representation (minor-symmetric; Case A)
///
/// * `m` pairs with the first two indices `ij`
/// * `n` equals `k`
/// * The symmetry is on `ij = ji`
/// 
/// ```text
/// Case A:
///      00    01    02
///    -----------------
/// 0 │ 00_0  00_1  00_2
/// 1 │ 11_0  11_1  11_2
/// 2 │ 22_0  22_1  22_2
/// 3 │ 01_0  01_1  01_2
/// 4 │ 12_0  12_1  12_2
/// 5 │ 02_0  02_1  02_2
///    -----------------
///      80    81    82
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::IJK_TO_MN_SYM_CASE_A;
/// assert_eq!(IJK_TO_MN_SYM_CASE_A[0][1][2], (3, 2));
/// ```
#[rustfmt::skip]
pub const IJK_TO_MN_SYM_CASE_A: [[[(usize, usize); 3]; 3]; 3] = [
    [
        [(0,0), (0,1), (0,2)], // 00_0, 00_1, 00_2
        [(3,0), (3,1), (3,2)], // 01_0, 01_1, 01_2
        [(5,0), (5,1), (5,2)], // 02_0, 02_1, 02_2
    ],
    [
        [(3,0), (3,1), (3,2)], // 01_0, 01_1, 01_2
        [(1,0), (1,1), (1,2)], // 11_0, 11_1, 11_2
        [(4,0), (4,1), (4,2)], // 12_0, 12_1, 12_2
    ],
    [
        [(5,0), (5,1), (5,2)], // 02_0, 02_1, 02_2
        [(4,0), (4,1), (4,2)], // 12_0, 12_1, 12_2
        [(2,0), (2,1), (2,2)], // 22_0, 22_1, 22_2
    ],
];

/// Maps (i,j,k) of Tensor3 to the (m,n)-th position in the matrix representation (minor-symmetric; Case B)
///
/// * `m` equals `i`
/// * `n` pairs with the last two indices `jk`
/// * The symmetry is on `jk = kj`
/// 
/// ```text
/// Case B:
///      00   01   02   03   04   05   06   07   08
///    ---------------------------------------------
/// 0 │ 0_00 0_11 0_22 0_01 0_12 0_02 0_01 0_12 0_02
/// 1 │ 1_00 1_11 1_22 1_01 1_12 1_02 1_01 1_12 1_02
/// 2 │ 2_00 2_11 2_22 2_01 2_12 2_02 2_01 2_12 2_02
///    ---------------------------------------------
///      20   21   22   23   24   25   26   27   28
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::IJK_TO_MN_SYM_CASE_B;
/// assert_eq!(IJK_TO_MN_SYM_CASE_B[0][1][2], (0, 4));
/// ```
#[rustfmt::skip]
pub const IJK_TO_MN_SYM_CASE_B: [[[(usize, usize); 3]; 3]; 3] = [
    [
        [(0,0),(0,3),(0,5)], // 0_00 0_01 0_02
        [(0,3),(0,1),(0,4)], // 0_01 0_11 0_12
        [(0,5),(0,4),(0,2)], // 0_02 0_12 0_22
    ],
    [
        [(1,0),(1,3),(1,5)], // 1_00 1_01 1_02
        [(1,3),(1,1),(1,4)], // 1_01 1_11 1_12
        [(1,5),(1,4),(1,2)], // 1_02 1_12 1_22
    ],
    [
        [(2,0),(2,3),(2,5)], // 2_00 2_01 2_02
        [(2,3),(2,1),(2,4)], // 2_01 2_11 2_12
        [(2,5),(2,4),(2,2)], // 2_02 2_12 2_22
    ],
];

//
// --- Fourth-order tensors ---
//

/// Maps the (m,n)-th position in the matrix representation to (i,j,k,l) of Tensor4
///
/// ```text
///      0  0   0  1   0  2    0  3   0  4   0  5    0  6   0  7   0  8
///    ----------------------------------------------------------------
/// 0 │ 00_00  00_11  00_22   00_01  00_12  00_02   00_10  00_21  00_20
/// 1 │ 11_00  11_11  11_22   11_01  11_12  11_02   11_10  11_21  11_20
/// 2 │ 22_00  22_11  22_22   22_01  22_12  22_02   22_10  22_21  22_20
///   │
/// 3 │ 01_00  01_11  01_22   01_01  01_12  01_02   01_10  01_21  01_20
/// 4 │ 12_00  12_11  12_22   12_01  12_12  12_02   12_10  12_21  12_20
/// 5 │ 02_00  02_11  02_22   02_01  02_12  02_02   02_10  02_21  02_20
///   │
/// 6 │ 10_00  10_11  10_22   10_01  10_12  10_02   10_10  10_21  10_20
/// 7 │ 21_00  21_11  21_22   21_01  21_12  21_02   21_10  21_21  21_20
/// 8 │ 20_00  20_11  20_22   20_01  20_12  20_02   20_10  20_21  20_20
///    ----------------------------------------------------------------
///      8  0   8  1   8  2    8  3   8  4   8  5    8  6   8  7   8  8
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::MN_TO_IJKL;
/// assert_eq!(MN_TO_IJKL[3][3], (0,1,0,1));
/// ```
#[rustfmt::skip]
pub const MN_TO_IJKL: [[(usize,usize,usize,usize); 9]; 9] = [
    [(0,0,0,0), (0,0,1,1), (0,0,2,2), (0,0,0,1), (0,0,1,2), (0,0,0,2), (0,0,1,0), (0,0,2,1), (0,0,2,0)], // 0
    [(1,1,0,0), (1,1,1,1), (1,1,2,2), (1,1,0,1), (1,1,1,2), (1,1,0,2), (1,1,1,0), (1,1,2,1), (1,1,2,0)], // 1
    [(2,2,0,0), (2,2,1,1), (2,2,2,2), (2,2,0,1), (2,2,1,2), (2,2,0,2), (2,2,1,0), (2,2,2,1), (2,2,2,0)], // 2
    [(0,1,0,0), (0,1,1,1), (0,1,2,2), (0,1,0,1), (0,1,1,2), (0,1,0,2), (0,1,1,0), (0,1,2,1), (0,1,2,0)], // 3
    [(1,2,0,0), (1,2,1,1), (1,2,2,2), (1,2,0,1), (1,2,1,2), (1,2,0,2), (1,2,1,0), (1,2,2,1), (1,2,2,0)], // 4
    [(0,2,0,0), (0,2,1,1), (0,2,2,2), (0,2,0,1), (0,2,1,2), (0,2,0,2), (0,2,1,0), (0,2,2,1), (0,2,2,0)], // 5
    [(1,0,0,0), (1,0,1,1), (1,0,2,2), (1,0,0,1), (1,0,1,2), (1,0,0,2), (1,0,1,0), (1,0,2,1), (1,0,2,0)], // 6
    [(2,1,0,0), (2,1,1,1), (2,1,2,2), (2,1,0,1), (2,1,1,2), (2,1,0,2), (2,1,1,0), (2,1,2,1), (2,1,2,0)], // 7
    [(2,0,0,0), (2,0,1,1), (2,0,2,2), (2,0,0,1), (2,0,1,2), (2,0,0,2), (2,0,1,0), (2,0,2,1), (2,0,2,0)], // 8
];

/// Maps (i,j,k,l) of Tensor4 to the (m,n)-th position in the matrix representation
///
/// ```text
///      0  0   0  1   0  2    0  3   0  4   0  5    0  6   0  7   0  8
///    ----------------------------------------------------------------
/// 0 │ 00_00  00_11  00_22   00_01  00_12  00_02   00_10  00_21  00_20
/// 1 │ 11_00  11_11  11_22   11_01  11_12  11_02   11_10  11_21  11_20
/// 2 │ 22_00  22_11  22_22   22_01  22_12  22_02   22_10  22_21  22_20
///   │
/// 3 │ 01_00  01_11  01_22   01_01  01_12  01_02   01_10  01_21  01_20
/// 4 │ 12_00  12_11  12_22   12_01  12_12  12_02   12_10  12_21  12_20
/// 5 │ 02_00  02_11  02_22   02_01  02_12  02_02   02_10  02_21  02_20
///   │
/// 6 │ 10_00  10_11  10_22   10_01  10_12  10_02   10_10  10_21  10_20
/// 7 │ 21_00  21_11  21_22   21_01  21_12  21_02   21_10  21_21  21_20
/// 8 │ 20_00  20_11  20_22   20_01  20_12  20_02   20_10  20_21  20_20
///    ----------------------------------------------------------------
///      8  0   8  1   8  2    8  3   8  4   8  5    8  6   8  7   8  8
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::IJKL_TO_MN;
/// assert_eq!(IJKL_TO_MN[0][1][0][1], (3,3));
/// ```
#[rustfmt::skip]
pub const IJKL_TO_MN: [[[[(usize, usize); 3]; 3]; 3]; 3] = [
    [
        [[(0,0), (0,3), (0,5)], [(0,6), (0,1), (0,4)], [(0,8), (0,7), (0,2)]], // 00[.][.]
        [[(3,0), (3,3), (3,5)], [(3,6), (3,1), (3,4)], [(3,8), (3,7), (3,2)]], // 01[.][.]
        [[(5,0), (5,3), (5,5)], [(5,6), (5,1), (5,4)], [(5,8), (5,7), (5,2)]], // 02[.][.]
    ],
    [
        [[(6,0), (6,3), (6,5)], [(6,6), (6,1), (6,4)], [(6,8), (6,7), (6,2)]], // 10[.][.]
        [[(1,0), (1,3), (1,5)], [(1,6), (1,1), (1,4)], [(1,8), (1,7), (1,2)]], // 11[.][.]
        [[(4,0), (4,3), (4,5)], [(4,6), (4,1), (4,4)], [(4,8), (4,7), (4,2)]], // 12[.][.]
    ],
    [
        [[(8,0), (8,3), (8,5)], [(8,6), (8,1), (8,4)], [(8,8), (8,7), (8,2)]], // 20[.][.]
        [[(7,0), (7,3), (7,5)], [(7,6), (7,1), (7,4)], [(7,8), (7,7), (7,2)]], // 21[.][.]
        [[(2,0), (2,3), (2,5)], [(2,6), (2,1), (2,4)], [(2,8), (2,7), (2,2)]], // 22[.][.]
    ],
];

/// Maps (i,j,k,l) of Tensor4 to the (m,n)-th position in the matrix representation (minor-symmetric version)
///
/// ```text
///      0  0   0  1   0  2    0  3   0  4   0  5
///    ------------------------------------------
/// 0 │ 00_00  00_11  00_22   00_01  00_12  00_02
/// 1 │ 11_00  11_11  11_22   11_01  11_12  11_02
/// 2 │ 22_00  22_11  22_22   22_01  22_12  22_02
///   │
/// 3 │ 01_00  01_11  01_22   01_01  01_12  01_02
/// 4 │ 12_00  12_11  12_22   12_01  12_12  12_02
/// 5 │ 02_00  02_11  02_22   02_01  02_12  02_02
///    ------------------------------------------
///      8  0   8  1   8  2    8  3   8  4   8  5
/// ```
///
/// # Examples
///
/// ```
/// use russell_tensor::IJKL_TO_MN_SYM;
/// assert_eq!(IJKL_TO_MN_SYM[0][1][0][1], (3,3));
/// ```
#[rustfmt::skip]
pub const IJKL_TO_MN_SYM: [[[[(usize, usize); 3]; 3]; 3]; 3] = [
    [
        [[(0,0), (0,3), (0,5)], [(0,3), (0,1), (0,4)], [(0,5), (0,4), (0,2)]], // 00[.][.]
        [[(3,0), (3,3), (3,5)], [(3,3), (3,1), (3,4)], [(3,5), (3,4), (3,2)]], // 01[.][.]
        [[(5,0), (5,3), (5,5)], [(5,3), (5,1), (5,4)], [(5,5), (5,4), (5,2)]], // 02[.][.]
    ],
    [
        [[(3,0), (3,3), (3,5)], [(3,3), (3,1), (3,4)], [(3,5), (3,4), (3,2)]], // 10[.][.]
        [[(1,0), (1,3), (1,5)], [(1,3), (1,1), (1,4)], [(1,5), (1,4), (1,2)]], // 11[.][.]
        [[(4,0), (4,3), (4,5)], [(4,3), (4,1), (4,4)], [(4,5), (4,4), (4,2)]], // 12[.][.]
    ],
    [
        [[(5,0), (5,3), (5,5)], [(5,3), (5,1), (5,4)], [(5,5), (5,4), (5,2)]], // 20[.][.]
        [[(4,0), (4,3), (4,5)], [(4,3), (4,1), (4,4)], [(4,5), (4,4), (4,2)]], // 21[.][.]
        [[(2,0), (2,3), (2,5)], [(2,3), (2,1), (2,4)], [(2,5), (2,4), (2,2)]], // 22[.][.]
    ],
];

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod tests {
    use super::{
        IJ_TO_M, IJ_TO_M_SYM, IJK_TO_MN_CASE_A, IJK_TO_MN_CASE_B, IJK_TO_MN_SYM_CASE_A, IJK_TO_MN_SYM_CASE_B,
        IJKL_TO_MN, IJKL_TO_MN_SYM, M_TO_IJ, MN_TO_IJK_CASE_A, MN_TO_IJK_CASE_B, MN_TO_IJKL, ONE_BY_3, SQRT_2,
        SQRT_2_BY_3, SQRT_3, SQRT_3_BY_2, SQRT_6, TWO_BY_3,
    };

    // #[test]
    // #[rustfmt::skip]
    // fn generate_mathematica_code() {
    //     print!("MNtoIJKL = {{");
    //     for m in 0..9 {
    //         if m > 0 { print!(","); } print!("{{");
    //         for n in 0..9 {
    //             if n > 0 { print!(","); }
    //             let (i, j, k, l) = MN_TO_IJKL[m][n];
    //             print!("{{{},{},{},{}}}", i + 1, j + 1, k + 1, l + 1);
    //         }
    //         print!("}}");
    //     }
    //     println!("}};");
    //     print!("IJKLtoMN = {{");
    //     for i in 0..3 {
    //         if i > 0 { print!(","); }
    //         print!("{{");
    //         for j in 0..3 {
    //             if j > 0 { print!(","); }
    //             print!("{{");
    //             for k in 0..3 {
    //                 if k > 0 { print!(","); } print!("{{");
    //                 for l in 0..3 {
    //                     if l > 0 {
    //                         print!(",");
    //                     }
    //                     let (m, n) = IJKL_TO_MN_SYM[i][j][k][l];
    //                     print!("{{{},{}}}", m + 1, n + 1);
    //                 }
    //                 print!("}}");
    //             }
    //             print!("}}");
    //         }
    //         print!("}}");
    //     }
    //     println!("}};");
    // }

    #[test]
    fn constants_are_correct() {
        assert_eq!(SQRT_2, 2_f64.sqrt());
        assert_eq!(SQRT_3, 3_f64.sqrt());
        assert_eq!(SQRT_6, 6_f64.sqrt());
        assert_eq!(SQRT_2_BY_3, (2_f64 / 3_f64).sqrt());
        assert_eq!(SQRT_3_BY_2, (3_f64 / 2_f64).sqrt());
        assert_eq!(ONE_BY_3, 1_f64 / 3_f64);
        assert_eq!(TWO_BY_3, 2_f64 / 3_f64);
    }

    #[test]
    fn tensor2_maps_are_correct() {
        //
        // Second-order tensors
        //
        // M_TO_IJ => IJ_TO_M and IJ_TO_M_SYM
        for m in 0..9 {
            let (i, j) = M_TO_IJ[m];
            assert_eq!(IJ_TO_M[i][j], m);
            let m_sym = match m {
                6 => 3,
                7 => 4,
                8 => 5,
                _ => m,
            };
            assert_eq!(IJ_TO_M_SYM[i][j], m_sym);
        }
    }

    #[test]
    fn tensor3_maps_are_correct_case_a() {
        //
        // Third-order tensors (Case A)
        //
        // MN_TO_IJK => IJK_TO_MN and IJK_TO_MN_SYM
        for m in 0..9 {
            let m_sym = match m {
                6 => 3,
                7 => 4,
                8 => 5,
                _ => m,
            };
            for n in 0..3 {
                let (i, j, k) = MN_TO_IJK_CASE_A[m][n];
                assert_eq!(IJK_TO_MN_CASE_A[i][j][k], (m, n));
                assert_eq!(IJK_TO_MN_SYM_CASE_A[i][j][k], (m_sym, n))
            }
        }
        // Compare with IJ_TO_M
        for i in 0..3 {
            for j in 0..3 {
                for k in 0..3 {
                    assert_eq!(IJK_TO_MN_CASE_A[i][j][k], (IJ_TO_M[i][j], k));
                }
            }
        }
    }

    #[test]
    fn tensor3_maps_are_correct_case_b() {
        //
        // Third-order tensors (Case B)
        //
        // MN_TO_IJK => IJK_TO_MN and IJK_TO_MN_SYM
        for m in 0..3 {
            for n in 0..9 {
                let (i, j, k) = MN_TO_IJK_CASE_B[m][n];
                assert_eq!(IJK_TO_MN_CASE_B[i][j][k], (m, n));
                let n_sym = match n {
                    6 => 3,
                    7 => 4,
                    8 => 5,
                    _ => n,
                };
                assert_eq!(IJK_TO_MN_SYM_CASE_B[i][j][k], (m, n_sym))
            }
        }
        // Compare with IJ_TO_M
        for i in 0..3 {
            for j in 0..3 {
                for k in 0..3 {
                    assert_eq!(IJK_TO_MN_CASE_B[i][j][k], (i, IJ_TO_M[j][k]));
                }
            }
        }
    }

    #[test]
    fn tensor4_maps_are_correct() {
        //
        // Fourth-order tensors
        //
        // MN_TO_IJKL => IJKL_TO_MN and IJKL_TO_MN_SYM
        for m in 0..9 {
            let m_sym = match m {
                6 => 3,
                7 => 4,
                8 => 5,
                _ => m,
            };
            for n in 0..9 {
                let (i, j, k, l) = MN_TO_IJKL[m][n];
                assert_eq!(IJKL_TO_MN[i][j][k][l], (m, n));
                let n_sym = match n {
                    6 => 3,
                    7 => 4,
                    8 => 5,
                    _ => n,
                };
                assert_eq!(IJKL_TO_MN_SYM[i][j][k][l], (m_sym, n_sym))
            }
        }
    }
}