1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

use crate::math::{DecomposedTensor, Matrix, Real, DIM};
use na::SMatrix;

#[cfg(not(feature = "std"))]
use na::ComplexField;

pub fn inv_exact(e: Real) -> Real {
    // We don't want to use any threshold here.
    if e == 0.0 {
        0.0
    } else {
        1.0 / e
    }
}

/// Computes the Lamé parameters (lambda, mu) from the young modulus and poisson ratio.
pub fn lame_lambda_mu(young_modulus: Real, poisson_ratio: Real) -> (Real, Real) {
    (
        young_modulus * poisson_ratio / ((1.0 + poisson_ratio) * (1.0 - 2.0 * poisson_ratio)),
        shear_modulus(young_modulus, poisson_ratio),
    )
}

pub fn shear_modulus(young_modulus: Real, poisson_ratio: Real) -> Real {
    young_modulus / (2.0 * (1.0 + poisson_ratio))
}

pub fn bulk_modulus(young_modulus: Real, poisson_ratio: Real) -> Real {
    young_modulus / (3.0 * (1.0 - 2.0 * poisson_ratio))
}

pub fn shear_modulus_from_lame(_lambda: Real, mu: Real) -> Real {
    mu
}

pub fn bulk_modulus_from_lame(lambda: Real, mu: Real) -> Real {
    lambda + 2.0 * mu / 3.0
}

pub fn solve_quadratic(a: Real, b: Real, c: Real) -> (Real, Real) {
    let discr_sqr = (b * b - 4.0 * a * c).sqrt();
    ((-b + discr_sqr) / (2.0 * a), (-b - discr_sqr) / (2.0 * a))
}

pub fn min_componentwise_quadratic_solve<const R: usize, const C: usize>(
    a: &SMatrix<Real, R, C>,
    b: &SMatrix<Real, R, C>,
    c: Real,
    sol_half_range: (Real, Real),
) -> Real {
    a.zip_map(b, |a, b| {
        let (mut s0, mut s1) = solve_quadratic(a, b, c);
        if s0 <= sol_half_range.0 {
            s0 = Real::MAX;
        }
        if s1 <= sol_half_range.0 {
            s1 = Real::MAX;
        }

        s0.min(s1)
    })
    .min()
    .min(sol_half_range.1)
}

pub fn spin_tensor(velocity_gradient: &Matrix<Real>) -> Matrix<Real> {
    (velocity_gradient - velocity_gradient.transpose()) * 0.5
}

pub fn strain_rate(velocity_gradient: &Matrix<Real>) -> Matrix<Real> {
    (velocity_gradient + velocity_gradient.transpose()) * 0.5
}

pub fn deviatoric_part(tensor: &Matrix<Real>) -> Matrix<Real> {
    DecomposedTensor::decompose(tensor).deviatoric_part
}
pub fn spherical_part(tensor: &Matrix<Real>) -> Real {
    tensor.trace() / (DIM as Real)
}