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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95

use crate::math::{Matrix, Real, Vector, DIM};

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

#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
#[derive(Copy, Clone, Debug, PartialEq)]
#[repr(C)]
pub struct DruckerPragerPlasticity {
    pub h0: Real,
    pub h1: Real,
    pub h2: Real,
    pub h3: Real,
    pub lambda: Real,
    pub mu: Real,
    pub only_active_when_failed: bool,
}

impl DruckerPragerPlasticity {
    pub fn new(young_modulus: Real, poisson_ration: Real) -> Self {
        let (lambda, mu) = crate::utils::lame_lambda_mu(young_modulus, poisson_ration);
        DruckerPragerPlasticity {
            h0: (35.0 as Real).to_radians(),
            h1: (9.0 as Real).to_radians(),
            h2: 0.2,
            h3: (10.0 as Real).to_radians(),
            lambda,
            mu,
            only_active_when_failed: false,
        }
    }

    pub fn project_deformation_gradient(
        &self,
        singular_values: Vector<Real>,
        log_vol_gain: Real,
        alpha: Real,
    ) -> Option<(Vector<Real>, Real)> {
        let d = DIM as Real;
        let strain = singular_values.map(|e| e.ln()) + Vector::repeat(log_vol_gain / d);
        let strain_trace = strain.sum();
        let deviatoric_strain = strain - Vector::repeat(strain_trace / d);

        if deviatoric_strain == Vector::zeros() || strain_trace > 0.0 {
            return Some((Vector::repeat(1.0), strain.norm()));
        }

        let gamma = deviatoric_strain.norm()
            + (d * self.lambda + 2.0 * self.mu) / (2.0 * self.mu) * strain_trace * alpha;
        if gamma <= 0.0 {
            return None;
        }

        let h = strain - gamma * deviatoric_strain.normalize();
        Some((h.map(|e| e.exp()), gamma))
    }

    fn alpha(&self, q: Real) -> Real {
        let angle = self.h0 + (self.h1 * q - self.h3) * (-self.h2 * q).exp();
        let s_angle = angle.sin();

        (2.0 as Real / 3.0).sqrt() * (2.0 * s_angle) / (3.0 - s_angle)
    }

    pub fn update_particle(
        &self,
        particle_phase: Real,
        particle_deformation_gradient: &mut Matrix<Real>,
        particle_plastic_deformation_gradient_det: &mut Real,
        particle_plastic_hardening: &mut Real,
        particle_log_vol_gain: &mut Real,
    ) {
        if self.only_active_when_failed && particle_phase != 0.0 {
            return;
        }

        let mut svd = particle_deformation_gradient.svd_unordered(true, true);

        let alpha = self.alpha(*particle_plastic_hardening);

        if let Some((new_singular_values, dq)) =
            self.project_deformation_gradient(svd.singular_values, *particle_log_vol_gain, alpha)
        {
            let prev_det = svd.singular_values.product();
            let new_det = new_singular_values.product();
            *particle_plastic_deformation_gradient_det *= prev_det / new_det;
            *particle_log_vol_gain += prev_det.ln() - new_det.ln();

            svd.singular_values = new_singular_values;
            *particle_plastic_hardening += dq; // TODO: Is it OK to store this here?
            *particle_deformation_gradient = svd.recompose().unwrap();
        } // Else, do nothing because we are inside of the yield surface.
    }
}