[wip] Derivatives

cpmech · Jul 2, 2023 · ac55bac · ac55bac
1 parent 1198ad7
commit ac55bac
Show file tree

Hide file tree

Showing 2 changed files with 108 additions and 181 deletions.
diff --git a/russell_tensor/src/order_high_derivatives.rs b/russell_tensor/src/order_high_derivatives.rs
@@ -1,5 +1,6 @@
 use crate::{
-    t2_dyad_t2, t2_odyad_t2, t2_qsd_t2, t2_ssd, Mandel, StrError, Tensor2, Tensor4, ONE_BY_3, SQRT_3, TOL_J2, TWO_BY_3,
+    t2_dyad_t2, t2_odyad_t2, t2_qsd_t2, t2_ssd, Deriv1InvariantJ2, Deriv1InvariantJ3, Mandel, StrError, Tensor2,
+    Tensor4, ONE_BY_3, SQRT_3, TOL_J2, TWO_BY_3,
 };
 use russell_lab::{mat_add, mat_mat_mul, mat_update};
 
@@ -322,6 +323,7 @@ impl Deriv2InvariantSigmaD {
         if jj2 > TOL_J2 {
             let a = 0.5 * SQRT_3 / f64::powf(jj2, 0.5);
             let b = 0.25 * SQRT_3 / f64::powf(jj2, 1.5);
+            Deriv1InvariantJ2::compute(&mut self.d1_jj2, sigma).unwrap();
             // sigma.deriv1_invariant_jj2(&mut self.d1_jj2).unwrap();
             deriv2_invariant_jj2(&mut self.d2_jj2, sigma).unwrap();
             t2_dyad_t2(d2, -b, &self.d1_jj2, &self.d1_jj2).unwrap();
@@ -337,6 +339,9 @@ pub struct Deriv2InvariantLode {
     /// auxiliary data to compute the second derivative of J3
     pub aux_jj3: Deriv2InvariantJ3,
 
+    /// TEMPORARY
+    pub temp: Deriv1InvariantJ3,
+
     /// auxiliary second-order tensor (Symmetric or Symmetric2D)
     pub tt: Tensor2,
 
@@ -370,6 +375,7 @@ impl Deriv2InvariantLode {
         }
         Ok(Deriv2InvariantLode {
             aux_jj3: Deriv2InvariantJ3::new(case).unwrap(),
+            temp: Deriv1InvariantJ3::new(case).unwrap(),
             tt: Tensor2::new(case),
             d1_jj2: Tensor2::new(case),
             d1_jj3: Tensor2::new(case),
@@ -421,6 +427,8 @@ impl Deriv2InvariantLode {
             let a = 1.5 * SQRT_3 / f64::powf(jj2, 1.5);
             let b = 2.25 * SQRT_3 / f64::powf(jj2, 2.5);
             let c = 5.625 * SQRT_3 / f64::powf(jj2, 3.5);
+            Deriv1InvariantJ2::compute(&mut self.d1_jj2, sigma).unwrap();
+            self.temp.compute(&mut self.d1_jj3, sigma).unwrap();
             // sigma.deriv1_invariant_jj2(&mut self.d1_jj2).unwrap();
             // sigma.deriv1_invariant_jj3(&mut self.d1_jj3, &mut self.tt).unwrap();
             deriv2_invariant_jj2(&mut self.d2_jj2, sigma).unwrap();
@@ -445,8 +453,8 @@ mod tests {
     use super::{Tensor2, Tensor4};
     use crate::{
         deriv2_invariant_jj2, deriv_inverse_tensor, deriv_inverse_tensor_sym, deriv_squared_tensor,
-        deriv_squared_tensor_sym, Deriv1InvariantJ2, Deriv2InvariantJ3, Deriv2InvariantLode, Deriv2InvariantSigmaD,
-        Mandel, SamplesTensor2, MN_TO_IJKL, SQRT_2,
+        deriv_squared_tensor_sym, Deriv1InvariantJ2, Deriv1InvariantJ3, Deriv1InvariantLode, Deriv1InvariantSigmaD,
+        Deriv2InvariantJ3, Deriv2InvariantLode, Deriv2InvariantSigmaD, Mandel, SamplesTensor2, MN_TO_IJKL, SQRT_2,
     };
     use russell_chk::{approx_eq, deriv_central5};
     use russell_lab::{mat_approx_eq, Matrix};
@@ -907,7 +915,7 @@ mod tests {
     struct ArgsNumDeriv2InvariantMandel {
         inv: Invariant, // option
         sigma: Tensor2, // temporary tensor
-        // d1: Tensor2,    // dInvariant/dσ (first derivative)
+        d1: Tensor2,    // dInvariant/dσ (first derivative)
         // s: Tensor2,     // deviator(σ)
         m: usize, // index of [dInvariant²/dσ⊗dσ]ₘₙ (Mandel components)
         n: usize, // index of [dInvariant²/dσ⊗dσ]ₘₙ (Mandel components)
@@ -916,36 +924,40 @@ mod tests {
     fn component_of_deriv1_inv_mandel(x: f64, args: &mut ArgsNumDeriv2InvariantMandel) -> f64 {
         let original = args.sigma.vec[args.n];
         args.sigma.vec[args.n] = x;
-        let res = match args.inv {
+        match args.inv {
             Invariant::J2 => {
-                let mut aux = Deriv1InvariantJ2::new(args.sigma.case()).unwrap();
-                aux.compute(&args.sigma).unwrap();
-                // args.d1.mirror(&aux.result).unwrap();
-                aux.result.vec[args.m]
-            } //args.sigma.deriv1_invariant_jj2(&mut args.d1).unwrap(),
-            Invariant::J3 => 0.0, //args.sigma.deriv1_invariant_jj3(&mut args.d1, &mut args.s).unwrap(),
+                Deriv1InvariantJ2::compute(&mut args.d1, &args.sigma).unwrap();
+                //args.sigma.deriv1_invariant_jj2(&mut args.d1).unwrap(),
+            }
+            Invariant::J3 => {
+                let mut aux = Deriv1InvariantJ3::new(args.sigma.case()).unwrap();
+                aux.compute(&mut args.d1, &args.sigma).unwrap();
+                //args.sigma.deriv1_invariant_jj3(&mut args.d1, &mut args.s).unwrap(),
+            }
             Invariant::SigmaD => {
+                Deriv1InvariantSigmaD::compute(&mut args.d1, &args.sigma)
+                    .unwrap()
+                    .unwrap();
                 // args.sigma.deriv1_invariant_sigma_d(&mut args.d1).unwrap().unwrap();
-                0.0
             }
             Invariant::Lode => {
+                let mut aux = Deriv1InvariantLode::new(args.sigma.case()).unwrap();
+                aux.compute(&mut args.d1, &args.sigma).unwrap().unwrap();
                 // args.sigma
                 //     .deriv1_invariant_lode(&mut args.d1, &mut args.s)
                 //     .unwrap()
                 //     .unwrap();
-                0.0
             }
         };
         args.sigma.vec[args.n] = original;
-        // args.d1.vec[args.m]
-        res
+        args.d1.vec[args.m]
     }
 
     fn numerical_deriv2_inv_sym_mandel(sigma: &Tensor2, inv: Invariant) -> Matrix {
         let mut args = ArgsNumDeriv2InvariantMandel {
             inv,
             sigma: Tensor2::new(Mandel::Symmetric),
-            // d1: Tensor2::new(Mandel::Symmetric),
+            d1: Tensor2::new(Mandel::Symmetric),
             // s: Tensor2::new(Mandel::Symmetric),
             m: 0,
             n: 0,