sparkl2d_core/utils/
physics.rs

1use crate::math::{DecomposedTensor, Matrix, Real, DIM};
2use na::SMatrix;
3
4#[cfg(not(feature = "std"))]
5use na::ComplexField;
6
7pub fn inv_exact(e: Real) -> Real {
8    // We don't want to use any threshold here.
9    if e == 0.0 {
10        0.0
11    } else {
12        1.0 / e
13    }
14}
15
16/// Computes the Lamé parameters (lambda, mu) from the young modulus and poisson ratio.
17/// https://encyclopediaofmath.org/wiki/Lam%C3%A9_constants
18pub fn lame_lambda_mu(young_modulus: Real, poisson_ratio: Real) -> (Real, Real) {
19    (
20        young_modulus * poisson_ratio / ((1.0 + poisson_ratio) * (1.0 - 2.0 * poisson_ratio)),
21        shear_modulus(young_modulus, poisson_ratio),
22    )
23}
24
25pub fn shear_modulus(young_modulus: Real, poisson_ratio: Real) -> Real {
26    young_modulus / (2.0 * (1.0 + poisson_ratio))
27}
28
29pub fn bulk_modulus(young_modulus: Real, poisson_ratio: Real) -> Real {
30    young_modulus / (3.0 * (1.0 - 2.0 * poisson_ratio))
31}
32
33pub fn shear_modulus_from_lame(_lambda: Real, mu: Real) -> Real {
34    mu
35}
36
37pub fn bulk_modulus_from_lame(lambda: Real, mu: Real) -> Real {
38    lambda + 2.0 * mu / 3.0
39}
40
41pub fn solve_quadratic(a: Real, b: Real, c: Real) -> (Real, Real) {
42    let discr_sqr = (b * b - 4.0 * a * c).sqrt();
43    ((-b + discr_sqr) / (2.0 * a), (-b - discr_sqr) / (2.0 * a))
44}
45
46pub fn min_componentwise_quadratic_solve<const R: usize, const C: usize>(
47    a: &SMatrix<Real, R, C>,
48    b: &SMatrix<Real, R, C>,
49    c: Real,
50    sol_half_range: (Real, Real),
51) -> Real {
52    a.zip_map(b, |a, b| {
53        let (mut s0, mut s1) = solve_quadratic(a, b, c);
54        if s0 <= sol_half_range.0 {
55            s0 = Real::MAX;
56        }
57        if s1 <= sol_half_range.0 {
58            s1 = Real::MAX;
59        }
60
61        s0.min(s1)
62    })
63    .min()
64    .min(sol_half_range.1)
65}
66
67pub fn spin_tensor(velocity_gradient: &Matrix<Real>) -> Matrix<Real> {
68    (velocity_gradient - velocity_gradient.transpose()) * 0.5
69}
70
71pub fn strain_rate(velocity_gradient: &Matrix<Real>) -> Matrix<Real> {
72    (velocity_gradient + velocity_gradient.transpose()) * 0.5
73}
74
75pub fn deviatoric_part(tensor: &Matrix<Real>) -> Matrix<Real> {
76    DecomposedTensor::decompose(tensor).deviatoric_part
77}
78pub fn spherical_part(tensor: &Matrix<Real>) -> Real {
79    tensor.trace() / (DIM as Real)
80}