Add eigenstrain (#70)

Author: dschwen

include/tensor_computes/FFTMechanics.h

@@ -65,6 +65,7 @@ class FFTMechanics : public TensorOperator<>
   const torch::Tensor * const _applied_macroscopic_strain;
 
   const bool _verbose;
+  const bool _accept_nonconverged;
 
   using TensorOperatorBase::_dim;
   using TensorOperatorBase::_domain;

include/tensor_computes/HyperElasticIsotropic.h

@@ -31,7 +31,11 @@ class HyperElasticIsotropic : public TensorOperator<>
   const torch::Tensor _tI4s;
   const torch::Tensor _tII;
 
+  /// deformation gradient tensor
   const torch::Tensor & _tF;
+  /// optional eigenstrain tensor
+  const torch::Tensor * _pFstar;
+
   const torch::Tensor & _tmu;
   const torch::Tensor & _tK;
   torch::Tensor & _tK4;

src/tensor_computes/FFTMechanics.C

@@ -42,6 +42,7 @@ FFTMechanics::validParams()
                                          "Applied macroscopic strain");
   params.addParam<TensorInputBufferName>("F", "F", "Deformation gradient tensor.");
   params.addParam<bool>("verbose", false, "Print non-linear residuals.");
+  params.addParam<bool>("accept_nonconverged", false, "Accept a non-converged solution");
   return params;
 }
 
@@ -69,7 +70,8 @@ FFTMechanics::FFTMechanics(const InputParameters & parameters)
     _applied_macroscopic_strain(isParamValid("applied_macroscopic_strain")
                                     ? &getInputBuffer("applied_macroscopic_strain")
                                     : nullptr),
-    _verbose(getParam<bool>("verbose"))
+    _verbose(getParam<bool>("verbose")),
+    _accept_nonconverged(getParam<bool>("accept_nonconverged"))
 {
   // Build projection tensor once
   const auto & q = _domain.getKGrid();
@@ -148,7 +150,7 @@ FFTMechanics::computeBuffer()
 
     // print nonlinear residual to the screen
     if (_verbose)
-      _console << "|R|=" << anorm << "\t|R/R0|=" << rnorm << '\n';
+      _console << "|R|=" << anorm << "\t|R/R0|=" << rnorm << std::endl;
 
     // check convergence
     if ((rnorm < _nl_rel_tol || anorm < _nl_abs_tol) && iiter > 0)
@@ -157,7 +159,17 @@ FFTMechanics::computeBuffer()
     iiter++;
 
     if (iiter > _nl_max_its)
-      paramError("nl_max_its",
-                 "Exceeded the maximum number of nonlinear iterations without converging.");
+    {
+      if (_accept_nonconverged)
+        paramWarning(
+            "nl_max_its", "Accepting non-converged solution |R|=", anorm, "\t|R/R0|=", rnorm);
+      else
+        paramError(
+            "nl_max_its",
+            "Exceeded the maximum number of nonlinear iterations without converging with |R|=",
+            anorm,
+            "\t|R/R0|=",
+            rnorm);
+    }
   }
 }

src/tensor_computes/HyperElasticIsotropic.C

@@ -17,6 +17,7 @@ HyperElasticIsotropic::validParams()
   InputParameters params = TensorOperator<>::validParams();
   params.addClassDescription("Hyperelastic isotropic constitutive model.");
   params.addRequiredParam<TensorInputBufferName>("F", "Deformation gradient tensor");
+  params.addParam<TensorInputBufferName>("Fstar", "Optional eigenstrain tensor");
   params.addRequiredParam<TensorInputBufferName>("mu", "Deformation gradient tensor");
   params.addRequiredParam<TensorInputBufferName>("K", "Deformation gradient tensor");
   params.addParam<TensorOutputBufferName>("tangent_operator", "dstressdstrain", "Stiffness tensor");
@@ -32,6 +33,7 @@ HyperElasticIsotropic::HyperElasticIsotropic(const InputParameters & parameters)
     _tI4s((_tI4 + _tI4rt) / 2.0),
     _tII(MooseTensor::dyad22(_tI, _tI)),
     _tF(getInputBuffer("F")),
+    _pFstar(isParamValid("Fstar") ? &getInputBuffer("Fstar") : nullptr),
     _tmu(getInputBuffer("mu")),
     _tK(getInputBuffer("K")),
     _tK4(getOutputBuffer("tangent_operator"))
@@ -43,10 +45,25 @@ HyperElasticIsotropic::computeBuffer()
 {
   using namespace MooseTensor;
 
+  const auto F_e = _pFstar ? dot22(_tF, torch::inverse(*_pFstar)) : _tF;
+
+  // Compute Green–Lagrange strain E = 0.5 (F_e^T * F_e - I)
+  const auto E = 0.5 * (dot22(trans2(F_e), F_e) - _tI);
+
+  // Build material stiffness tensor
   const auto C4 = _tK.reshape(_domain.getValueShape({1, 1, 1, 1})) * _tII +
                   2. * _tmu.reshape(_domain.getValueShape({1, 1, 1, 1})) * (_tI4s - 1. / 3. * _tII);
-  const auto S = ddot42(C4, .5 * (dot22(trans2(_tF), _tF) - _tI));
 
-  _u = dot22(_tF, S);
-  _tK4 = dot24(S, _tI4) + ddot44(ddot44(_tI4rt, dot42(dot24(_tF, C4), trans2(_tF))), _tI4rt);
+  // Second Piola-Kirchhoff stress: S = C : E
+  const auto S = ddot42(C4, E);
+
+  // First Piola-Kirchhoff stress: P = F_e * S
+  _u = dot22(F_e, S);
+
+  _tK4 = dot24(S, _tI4) + ddot44(ddot44(_tI4rt, dot42(dot24(F_e, C4), trans2(F_e))), _tI4rt);
+
+  // Consistent tangent: K_ijkl = F^e_{im} C_{jlmn} F^e_{kn} + delta_{ik} S_{jl}
+  // const auto K_geo = torch::einsum("...im,...jlmn,...kn->...ijkl", {F_e, C4, F_e});
+  // const auto K_stress = torch::einsum("ik,jl,...jl->...ijkl", {_ti, _ti, S});
+  // _tK4 = K_geo + K_stress;
 }