sparkl2d_core/utils/
physics.rs1use 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 if e == 0.0 {
10 0.0
11 } else {
12 1.0 / e
13 }
14}
15
16pub 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}