#123 using plastic Jacobian_mult

+ fixed error in over computation for MCC

src/materials/RedbackMechMaterial.C

@@ -681,10 +681,9 @@ RedbackMechMaterial::returnMap(const RankTwoTensor & sig_old,
   const Real tol1 = 1e-10; // TODO: expose to user interface and/or make the tolerance relative
   const Real tol3 = 1e-6;  // TODO: expose to user interface and/or make the tolerance relative
   Real err3 = 1.1 * tol3;
-  bool is_plastic; // is this point in plastic regime or not?
-  bool is_first_plastic_determined = false;
-  bool is_first_plastic; // is_plastic the first time it's called
-  Real s;                // curvilinear arc-length between (p,q) and (p_y,q_y)
+  bool is_plastic = false; // is this point in plastic regime or not?
+  Real s = 0; // curvilinear arc-length between (p,q) and (p_y,q_y)
+  RankFourTensor dr_dsig_inv = E_ijkl;
 
   Real eqvpstrain = std::pow(2.0 / 3.0, 0.5) * dp.L2norm();
   Real yield_stress = getYieldStress(eqvpstrain);
@@ -726,12 +725,8 @@ RedbackMechMaterial::returnMap(const RankTwoTensor & sig_old,
     Real p = sig_new.trace() / 3.0;
     Real q = getSigEqv(sig_new);
     get_py_qy_damaged(p, q, p_y, q_y, yield_stress, is_plastic, s);
-    if (!is_first_plastic_determined)
-    {
-      is_first_plastic = is_plastic;
-      is_first_plastic_determined = true;
-    }
-    if (is_first_plastic)
+
+    if (is_plastic)
     {
       Real flow_incr = getFlowIncrement(q, p, q_y, p_y, yield_stress, s);
 
@@ -749,7 +744,7 @@ RedbackMechMaterial::returnMap(const RankTwoTensor & sig_old,
         // Jacobian = d(residual)/d(sigma)
         RankFourTensor dr_dsig;
         getJac(sig_new, E_ijkl, flow_incr, q, p, p_y, q_y, yield_stress, s, dr_dsig);
-        RankFourTensor dr_dsig_inv = dr_dsig.invSymm();
+        dr_dsig_inv = dr_dsig.invSymm();
         RankTwoTensor ddsig = -dr_dsig_inv * resid; // Newton Raphson
         delta_dp -= E_ijkl.invSymm() * ddsig;       // Update increment of plastic rate of deformation tensor
         sig_new += ddsig;                           // Update stress
@@ -788,6 +783,9 @@ RedbackMechMaterial::returnMap(const RankTwoTensor & sig_old,
     throw MooseException("Constitutive Error-Too many iterations in Hardness "
                          "Update:Reduce time increment.\n"); // Convergence failure
 
+  if (is_plastic && _t_step > 1) // keep elasticity jacobian_mult for very first step
+    _Jacobian_mult[ _qp ] = dr_dsig_inv;
+
   dp = dpn; // Plastic rate of deformation tensor in unrotated configuration
   sig = sig_new;
 }

src/utils/Ellipse.C

@@ -176,7 +176,7 @@ Ellipse::getYieldPointCC(
   Real const m, Real const p_c, Real const y0, Real const y1, Real & x0, Real & x1, Real & s, Real shift)
 {
   Real t; // curvilinear "time"
-  Real alpha, beta, gamma, norm_n;
+  Real alpha, beta, gamma;
   if (y1 < 0)
   {
     // algorithm below only works for y1 > 0, so use symmetry
@@ -244,20 +244,16 @@ Ellipse::getYieldPointCC(
   beta = 2 * (1 - 1 / std::pow(m, 2));
   // If all arguments |z|<1, we can use hypergeometric function 2F1 (not accounting for shift!)
   // (as we don't have a good implementation for |z|>1)
-  norm_n = std::sqrt(std::pow(2 * y0 - p_c - 2 * shift, 2) / 3 + 6 * std::pow(y1, 2) / std::pow(m, 4));
-  if (std::abs(alpha) < 1 && std::abs(alpha * std::exp(beta * 2 * t / norm_n)) < 1)
+  if (std::abs(alpha) < 1 && std::abs(alpha * std::exp(beta * 2 * t)) < 1)
   {
     s = std::fabs(gamma) / (beta - 2.0) *
         (std::exp(-2 * t) *
-           (2 * std::sqrt(1 + alpha * std::exp(beta * 2 * t / norm_n)) -
-            beta * hyp2f1(0.5, -1 / beta, (beta - 1.0) / beta, -alpha * std::exp(beta * 2 * t / norm_n))) -
+           (2 * std::sqrt(1 + alpha * std::exp(beta * 2 * t)) -
+            beta * hyp2f1(0.5, -1 / beta, (beta - 1.0) / beta, -alpha * std::exp(beta * 2 * t))) -
          2 * std::sqrt(1 + alpha) + beta * hyp2f1(0.5, -1 / beta, (beta - 1.0) / beta, -alpha));
   }
   else
   {
-    // Real test = alpha*std::exp(beta*2*t/norm_n);
-    // std::cout << "getYieldPointCC using sum of segments with shift=" << shift << ", alpha=" << alpha << ",
-    // alpha*std::exp(beta*2*t/norm_n)="<< test << std::endl;
     int n_iter = 100; // TODO: is this value good enough even very far from the ellipse?
     Real t_old = 0;
     Real t_new = 0;

unit/src/EllipseTest.C

@@ -831,12 +831,12 @@ EllipseTest::distanceCCTestMajorAxisVertical()
   // top right quadrant
   d = Ellipse::distanceCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-0.5, /*y1=*/2.0, x2, y2, 0.);
   Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-0.5, /*y1=*/2.0, x3, y3, s, 0.);
-  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.2094051, d, 1e-5);
-  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.209475269026, s, 1e-5);
-  CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.650317, x2, 1e-5);
-  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.854209, y2, 1e-5);
-  CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.653371450454, x3, 1e-10);
-  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.85734897716, y3, 1e-10);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE("major axis vertical, wrong d", 0.2094051, d, 1e-5);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE("major axis vertical, wrong s", 0.209475269026, s, 1e-5);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE("major axis vertical, wrong x2", -0.650317, x2, 1e-5);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE("major axis vertical, wrong y2", 1.854209, y2, 1e-5);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE("major axis vertical, wrong x3", -0.653371450454, x3, 1e-10);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE("major axis vertical, wrong y3", 1.85734897716, y3, 1e-10);
   // top left quadrant
   d = Ellipse::distanceCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-2.5, /*y1=*/2.0, x2, y2, 0.);
   Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-2.5, /*y1=*/2.0, x3, y3, s, 0.);