Elimnated use of stiffness tensor.

MaterialLib/SolidModels/MFront/Lubby2.mfront

@@ -18,7 +18,6 @@ rheological model
 @Parameter local_zero_tolerance = 1.e-14;
 
 @ModellingHypotheses{".+"};
-@RequireStiffnessTensor<UnAltered>;
 
 // Intercept of yield function
 @MaterialProperty real GK0;
@@ -40,14 +39,21 @@ epsK.setEntryName("KelvinStrain");
 @StateVariable Stensor epsM;
 epsM.setEntryName("MaxwellStrain");
 
+//! Second Lamé coefficient
+@LocalVariable stress mu;
+
+@InitLocalVariables
+{
+    mu = computeMu(young, nu);
+}
+
 @Integrator
 {
     // Implicit system
     constexpr auto id4 = Stensor4::Id();
     constexpr auto Pdev = Stensor4::K();
 
-    const auto sigma_eff =
-        std::max(sigmaeq(sig), local_zero_tolerance * D(0, 0));
+    const auto sigma_eff = std::max(sigmaeq(sig), local_zero_tolerance * mu);
 
     const auto s = deviator(sig);
     const auto k = deviator(epsK + depsK);
@@ -66,8 +72,11 @@ epsM.setEntryName("MaxwellStrain");
     // calculate Maxwell strain residual
     fepsM = depsM - dt / (2.0 * etaM) * s;
 
+    // Pdev * D
+    const auto Pdev_D = 2 * mu * Pdev;
+
     // Jacobian
-    const auto dsigma_eff_ddeel = 3. / (2. * sigma_eff) * s | (Pdev * D);
+    const auto dsigma_eff_ddeel = 3. / (2. * sigma_eff) * s | (Pdev_D);
     const auto detaK_ddeel = etaK * mvK * dsigma_eff_ddeel;
     const auto detaM_ddeel = etaM * mvM * dsigma_eff_ddeel;
     const auto dGK_ddeel = GK * mK * dsigma_eff_ddeel;
@@ -77,12 +86,12 @@ epsM.setEntryName("MaxwellStrain");
     dfeel_ddepsM = id4;
 
     dfepsK_ddeel = (dt / (2.0 * etaK * etaK) * s_m_2G_eK ^ detaK_ddeel) -
-                   (dt / (2.0 * etaK) * Pdev * D) + (dt / etaK * k ^ dGK_ddeel);
+                   (dt / (2.0 * etaK) * Pdev_D) + (dt / etaK * k ^ dGK_ddeel);
     dfepsK_ddepsK = (1. + dt * GK / etaK) * id4;
     // dfepsK_ddepsM is all zero --> nothing to do
 
     dfepsM_ddeel = (dt / (2.0 * etaM * etaM) * s ^ detaM_ddeel) -
-                   (dt / (2.0 * etaM) * Pdev * D);
+                   (dt / (2.0 * etaM) * Pdev_D);
     // dfepsM_ddepsK is all zero --> nothing to do
     dfepsM_ddepsM = id4;
 }

MaterialLib/SolidModels/MFront/Lubby2mod.mfront

@@ -2,7 +2,8 @@
 @DSL Implicit;
 @Behaviour Lubby2mod;
 @Author Thomas Nagel, Thomas Helfer;
-@Description {Lubby2 model for stationary
+@Description
+{Lubby2 model for stationary
 and transient creep based on Burgers
 rheological model.
 Reduced implementation.
@@ -16,7 +17,6 @@ Reduced implementation.
 
 @ModellingHypotheses{".+"};
 @Brick StandardElasticity;
-@RequireStiffnessTensor<UnAltered>;
 
 // Intercept of yield function
 @MaterialProperty real GK0;
@@ -40,41 +40,53 @@ epsK.setEntryName("KelvinStrain");
 //! increment of the Maxwell strain
 @LocalVariable Stensor depsM;
 
-@Integrator {
-
-  const auto seq = sigmaeq(sig);
-  const auto iseq = 1. / std::max(seq, 1.e-14 * D(0, 0));
-  const auto s = deviator(sig);
-  constexpr auto Pdev = Stensor4::K();
-
-  const auto etaK = etaK0 * std::exp(mvK * seq);
-  const auto etaM = etaM0 * std::exp(mvM * seq);
-  const auto GK = GK0 * std::exp(mK * seq);
-
-  // calculate Kelvin strain residual
-  const auto etaK_impl = etaK + dt * GK * theta;
-  depsK = dt / (2. * etaK_impl) * (s - 2 * GK * epsK);
-
-  // calculate Maxwell strain residual
-  depsM = dt / (2 * etaM) * s;
-
-  // residuals
-  feel += depsM + depsK;
-
-  // Jacobian
-  const auto dseq_ddeel = 3. * iseq / 2. * s | (Pdev * D);
-  const auto detaK_ddeel = etaK * mvK * dseq_ddeel;
-  const auto detaM_ddeel = etaM * mvM * dseq_ddeel;
-  const auto dGK_ddeel = GK * mK * dseq_ddeel;
-  //
-  dfeel_ddeel +=
-      -(dt / (2 * etaK_impl * etaK_impl) * (s - 2 * GK * epsK) ^ detaK_ddeel) +
-      (dt / (2 * etaK_impl) * Pdev * D) - dt / etaK_impl * (epsK ^ dGK_ddeel) -
-      dt / (2 * etaK_impl * etaK_impl) * dt * theta *
-          ((s - 2 * GK * epsK) ^ dGK_ddeel) -
-      (dt / (2 * etaM * etaM) * s ^ detaM_ddeel) + (dt / (2 * etaM) * Pdev * D);
+//! Second Lamé coefficient
+@LocalVariable stress mu;
+
+@InitLocalVariables
+{
+    mu = computeMu(young, nu);
+}
+
+@Integrator
+{
+    const auto seq = sigmaeq(sig);
+    const auto iseq = 1. / std::max(seq, 1.e-14 * mu);
+    const auto s = deviator(sig);
+    constexpr auto Pdev = Stensor4::K();
+
+    const auto etaK = etaK0 * std::exp(mvK * seq);
+    const auto etaM = etaM0 * std::exp(mvM * seq);
+    const auto GK = GK0 * std::exp(mK * seq);
+
+    // calculate Kelvin strain residual
+    const auto etaK_impl = etaK + dt * GK * theta;
+    depsK = dt / (2. * etaK_impl) * (s - 2 * GK * epsK);
+
+    // calculate Maxwell strain residual
+    depsM = dt / (2 * etaM) * s;
+
+    // residuals
+    feel += depsM + depsK;
+
+    // Jacobian
+    // Pdev_D
+    const auto Pdev_D = 2 * mu * Pdev;
+    const auto dseq_ddeel = 3. * iseq / 2. * s | (Pdev_D);
+    const auto detaK_ddeel = etaK * mvK * dseq_ddeel;
+    const auto detaM_ddeel = etaM * mvM * dseq_ddeel;
+    const auto dGK_ddeel = GK * mK * dseq_ddeel;
+    //
+    dfeel_ddeel +=
+        -(dt / (2 * etaK_impl * etaK_impl) * (s - 2 * GK * epsK) ^
+          detaK_ddeel) +
+        (dt / (2 * etaK_impl) * Pdev_D) - dt / etaK_impl * (epsK ^ dGK_ddeel) -
+        dt / (2 * etaK_impl * etaK_impl) * dt * theta *
+            ((s - 2 * GK * epsK) ^ dGK_ddeel) -
+        (dt / (2 * etaM * etaM) * s ^ detaM_ddeel) + (dt / (2 * etaM) * Pdev_D);
 }
 
-@UpdateAuxiliaryStateVariables {
-  epsK += depsK;
-} // end of @UpdateAuxiliaryStateVariables
+@UpdateAuxiliaryStateVariables
+{
+    epsK += depsK;
+}  // end of @UpdateAuxiliaryStateVariables