From c197722ed267bafbb3bf45dfc25f6d9b004145ba Mon Sep 17 00:00:00 2001
From: Prithivirajan Veerappan <veerappan.prithivirajan@inl.gov>
Date: Tue, 30 Jul 2024 16:31:15 -0600
Subject: [PATCH] Addressed review comments - Minor changes ref #184

---
 .../materials/DamagePlasticityStressUpdate.h  |  1 -
 src/materials/DamagePlasticityStressUpdate.C  | 22 ++++++++++++++-----
 2 files changed, 17 insertions(+), 6 deletions(-)

diff --git a/include/materials/DamagePlasticityStressUpdate.h b/include/materials/DamagePlasticityStressUpdate.h
index c7cf6b6c8..ee4b66810 100644
--- a/include/materials/DamagePlasticityStressUpdate.h
+++ b/include/materials/DamagePlasticityStressUpdate.h
@@ -231,7 +231,6 @@ class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate
    * @return dr_dstress (dflowpotential_dstress)
    */
   virtual void dflowPotential_dstress(const std::vector<Real> & stress_params,
-                                      const std::vector<Real> & intnl,
                                       std::vector<std::vector<Real>> & dr_dstress) const;
   /**
    * This function calculates the derivative of the flow potential with the damage states
diff --git a/src/materials/DamagePlasticityStressUpdate.C b/src/materials/DamagePlasticityStressUpdate.C
index 0cfc8883f..969830bb2 100644
--- a/src/materials/DamagePlasticityStressUpdate.C
+++ b/src/materials/DamagePlasticityStressUpdate.C
@@ -130,6 +130,16 @@ DamagePlasticityStressUpdate::initQpStatefulProperties()
       used, the following commented lines show several different options. Some other options are
       still being considered. In this code, we define the element length as the cube root of the
       element volume */
+
+  // if (_current_elem->n_vertices() < 3)
+  //   _ele_len[_qp] = _current_elem->length(0, 1);
+  // else if (_current_elem->n_vertices() < 5)
+  //   _ele_len[_qp] = (_current_elem->length(0, 1) + _current_elem->length(1, 2)) / 2.;
+  // else
+  //   _ele_len[_qp] =
+  //       (_current_elem->length(0, 1) + _current_elem->length(1, 2) + _current_elem->length(0, 4))
+  //       / 3.;
+
   _ele_len[_qp] = std::cbrt(_current_elem->volume());
 
   _gt[_qp] = _FEt / _ele_len[_qp];
@@ -315,9 +325,7 @@ DamagePlasticityStressUpdate::flowPotential(const std::vector<Real> & stress_par
 
 void
 DamagePlasticityStressUpdate::dflowPotential_dstress(
-    const std::vector<Real> & stress_params,
-    const std::vector<Real> & /* intnl */,
-    std::vector<std::vector<Real>> & dr_dstress) const
+    const std::vector<Real> & stress_params, std::vector<std::vector<Real>> & dr_dstress) const
 {
   Real J2 = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0)
                 .secondInvariant();
@@ -522,7 +530,9 @@ DamagePlasticityStressUpdate::fbar(const Real & f0,
   Real v = sqrt_phi;
   Real u = (1 + a) / a - sqrt_phi / a;
   Real cal_fbar = f0 * std::pow(u, exponent) * v;
-  return (u > 0) ? cal_fbar : 1.E-6;
+  return (u > 0) ? cal_fbar : 1.E-6; // The minimum value for the strength parameter is 1.E-6, as a
+                                     // zero value can cause numerical instabilities due to
+                                     // derivatives and use of logarithmic functions.
 }
 
 Real
@@ -552,7 +562,9 @@ DamagePlasticityStressUpdate::f(const Real & f0, const Real & a, const Real & ka
   Real v = phi;
   Real u = (1 + a) * sqrt_phi;
   Real cal_f = (f0 / a) * (u - v);
-  return (u > v) ? cal_f : 1.E-6;
+  return (u > v) ? cal_f : 1.E-6; // The minimum value for the strength parameter is 1.E-6, as a
+                                  // zero value can cause numerical instabilities due to derivatives
+                                  // and use of logarithmic functions.
 }
 
 Real