Inelastic heat can be used in a MultiApp as well

include/kernels/LynxHeatSources.h

@@ -46,7 +46,7 @@ protected:
   const MaterialProperty<Real> & _radiogenic_heat;
   const MaterialProperty<Real> & _rho_b;
   const MaterialProperty<Real> & _dinvrho_dtemp;
-  const MaterialProperty<Real> & _inelastic_heat;
+  const MaterialProperty<Real> & _inelastic_heat_mat;
   const MaterialProperty<Real> & _adiabatic_heat;
   const MaterialProperty<Real> & _damage_heat;
   // const MaterialProperty<Real> & _dinelastic_heat_dtemp;
@@ -54,6 +54,8 @@ protected:
   const bool _coupled_disp;
   unsigned int _ndisp;
   std::vector<unsigned int> _disp_var;
+  const bool _coupled_inelastic_heat;
+  const VariableValue & _inelastic_heat;
 };
 
 #endif // LYNXHEATSOURCES_H

include/kernels/LynxMass.h

@@ -45,12 +45,12 @@ protected:
   virtual Real computeQpJacobian();
   virtual Real computeQpOffDiagJacobian(unsigned jvar);
 
+  const Real _penalty;
   PenaltyType _penalty_type;
 
   unsigned _ndisp;
   std::vector<unsigned> _disp_var;
 
-  const MaterialProperty<Real> & _penalty;
   const MaterialProperty<RankTwoTensor> & _strain_increment;
 };
 

include/kernels/LynxSolidMomentum.h

@@ -39,6 +39,7 @@ public:
   virtual void computeResidual() override;
   virtual void computeJacobian() override;
   virtual void computeOffDiagJacobian(MooseVariableFEBase & jvar) override;
+  using Kernel::computeOffDiagJacobian;
 
 protected:
   virtual Real computeQpResidual() override;

src/kernels/LynxHeatSources.C

@@ -35,6 +35,7 @@ validParams<LynxHeatSources>()
       "coeff_shear_heating", 0.0, "The coefficient in front of the shear heating generation.");
   params.addCoupledVar("displacements",
                        "The string of displacements suitable for the problem statement");
+  params.addCoupledVar("inelastic_heat", "The auxiliary variable holding the inelastic heat value for running in a subApp.");
   return params;
 }
 
@@ -45,13 +46,15 @@ LynxHeatSources::LynxHeatSources(const InputParameters & parameters)
     _radiogenic_heat(getDefaultMaterialProperty<Real>("radiogenic_heat_production")),
     _rho_b(getDefaultMaterialProperty<Real>("bulk_density")),
     _dinvrho_dtemp(getDefaultMaterialProperty<Real>("dinvrho_dtemp")),
-    _inelastic_heat(getDefaultMaterialProperty<Real>("inelastic_heat")),
+    _inelastic_heat_mat(getDefaultMaterialProperty<Real>("inelastic_heat")),
     _adiabatic_heat(getDefaultMaterialProperty<Real>("adiabatic_heat")),
     _damage_heat(getDefaultMaterialProperty<Real>("damage_heat")),
     _dinelastic_heat_dstrain(getDefaultMaterialProperty<RankTwoTensor>("dinelastic_heat_dstrain")),
     _coupled_disp(isCoupled("displacements")),
     _ndisp(_coupled_disp ? coupledComponents("displacements") : 0),
-    _disp_var(_ndisp)
+    _disp_var(_ndisp),
+    _coupled_inelastic_heat(isCoupled("inelastic_heat")),
+    _inelastic_heat(_coupled_inelastic_heat ? coupledValue("inelastic_heat") : _zero)
 {
   for (unsigned int i = 0; i < _ndisp; ++i)
     _disp_var[i] = coupled("displacements", i);
@@ -70,7 +73,11 @@ Real
 LynxHeatSources::computeQpResidual()
 {
   Real Hr = _radiogenic_heat[_qp];
-  Real Hs = _coeff_Hs * _inelastic_heat[_qp];
+  Real Hs = _coeff_Hs;
+  if (_coupled_inelastic_heat)
+    Hs *= _inelastic_heat[_qp];
+  else
+    Hs *= _inelastic_heat_mat[_qp];
   Real Ha = _adiabatic_heat[_qp];
   Real Hd = _damage_heat[_qp];
   Real heat_sources = Hr + Hs + Ha + Hd;
@@ -88,7 +95,11 @@ Real
 LynxHeatSources::computeQpJacobian()
 {
   Real Hr = _radiogenic_heat[_qp];
-  Real Hs = _coeff_Hs * _inelastic_heat[_qp];
+  Real Hs = _coeff_Hs;
+  if (_coupled_inelastic_heat)
+    Hs *= _inelastic_heat[_qp];
+  else
+    Hs *= _inelastic_heat_mat[_qp];
   Real Ha = _adiabatic_heat[_qp];
   Real Hd = _damage_heat[_qp];
   Real heat_sources = Hr + Hs + Ha + Hd;

src/kernels/LynxMass.C

@@ -29,6 +29,7 @@ validParams<LynxMass>()
   params.addClassDescription("Divergence of solid velocity for incompressible Stoke flow.");
   params.addRequiredCoupledVar(
       "displacements", "The string of displacements variables suitable for the problem statement.");
+  params.addParam<Real>("penalty", 0.0, "The value of the penalty.");
   MooseEnum penalty_type_options("linear=0 laplace=1", "linear");
   params.addParam<MooseEnum>(
       "penalty_type", penalty_type_options, "The type of penalty formulation.");
@@ -38,10 +39,10 @@ validParams<LynxMass>()
 
 LynxMass::LynxMass(const InputParameters & parameters)
   : DerivativeMaterialInterface<Kernel>(parameters),
+    _penalty(getParam<Real>("penalty")),
     _penalty_type((PenaltyType)(int)parameters.get<MooseEnum>("penalty_type")),
     _ndisp(coupledComponents("displacements")),
     _disp_var(_ndisp),
-    _penalty(getDefaultMaterialProperty<Real>("penalty")),
     _strain_increment(getDefaultMaterialProperty<RankTwoTensor>("strain_increment"))
 {
   for (unsigned int i = 0; i < _ndisp; ++i)
@@ -55,7 +56,7 @@ LynxMass::LynxMass(const InputParameters & parameters)
 Real
 LynxMass::computeQpResidual()
 {
-  Real one_on_penalty = (_penalty[_qp] != 0) ? 1.0 / _penalty[_qp] : 0.0;
+  Real one_on_penalty = (_penalty != 0) ? 1.0 / _penalty : 0.0;
   Real res = -_strain_increment[_qp].trace() * _test[_i][_qp];
   switch (_penalty_type)
   {
@@ -77,7 +78,7 @@ LynxMass::computeQpResidual()
 Real
 LynxMass::computeQpJacobian()
 {
-  Real one_on_penalty = (_penalty[_qp] != 0) ? 1.0 / _penalty[_qp] : 0.0;
+  Real one_on_penalty = (_penalty != 0) ? 1.0 / _penalty : 0.0;
   Real jac = 0.0;
 
   switch (_penalty_type)

src/materials/LynxDeformationBase.C

@@ -675,10 +675,11 @@ LynxDeformationBase::updateSpinTangentModulus()
   RankTwoTensor strain_el_dev_old = _elastic_strain_old[_qp].deviatoric();
   RankTwoTensor Id = RankTwoTensor(RankTwoTensor::initIdentity);
 
-  _tangent_modulus[_qp] +=
-      _G[_qp] * (Id.mixedProductIkJl(strain_el_dev_old.transpose()) -
-                 Id.mixedProductIlJk(strain_el_dev_old.transpose())) -
-      _G[_qp] * (strain_el_dev_old.mixedProductIkJl(Id) - strain_el_dev_old.mixedProductIlJk(Id));
+  if (_G[_qp] != 0.0)
+    _tangent_modulus[_qp] +=
+        _G[_qp] * (Id.mixedProductIkJl(strain_el_dev_old.transpose()) -
+                   Id.mixedProductIlJk(strain_el_dev_old.transpose())) -
+        _G[_qp] * (strain_el_dev_old.mixedProductIkJl(Id) - strain_el_dev_old.mixedProductIlJk(Id));
 }
 
 void