diff --git a/MaterialLib/SolidModels/Lubby2-impl.h b/MaterialLib/SolidModels/Lubby2-impl.h
index a52ace4bdcd956d8a6e3519fbf143b1fac0f66e3..6f83f05140e95d51d4bf7199338189e4fd12074c 100644
--- a/MaterialLib/SolidModels/Lubby2-impl.h
+++ b/MaterialLib/SolidModels/Lubby2-impl.h
@@ -34,6 +34,11 @@ bool Lubby2<DisplacementDim>::computeConstitutiveRelation(
nullptr);
MaterialStateVariables& state =
static_cast<MaterialStateVariables&>(material_state_variables);
+ state.loadState();
+
+ auto local_lubby2_properties =
+ detail::LocalLubby2Properties<DisplacementDim>{t, x, _mp};
+
// calculation of deviatoric parts
auto const& P_dev = Invariants::deviatoric_projection;
KelvinVector const epsd_i = P_dev * eps;
@@ -43,7 +48,7 @@ bool Lubby2<DisplacementDim>::computeConstitutiveRelation(
// Calculate effective stress and update material properties
double sig_eff = Invariants::equivalentStress(sigd_j);
- updateBurgersProperties(t, x, sig_eff * state.GM, state);
+ local_lubby2_properties.update(sig_eff);
using LocalJacobianMatrix =
Eigen::Matrix<double, KelvinVectorSize * 3, KelvinVectorSize * 3,
@@ -71,15 +76,15 @@ bool Lubby2<DisplacementDim>::computeConstitutiveRelation(
LocalJacobianMatrix K_loc;
auto const update_residual = [&](LocalResidualVector& residual) {
- calculateResidualBurgers(dt, epsd_i, sigd_j, state.eps_K_j,
- state.eps_K_t, state.eps_M_j,
- state.eps_M_t, residual, state);
+ calculateResidualBurgers(
+ dt, epsd_i, sigd_j, state.eps_K_j, state.eps_K_t, state.eps_M_j,
+ state.eps_M_t, residual, local_lubby2_properties);
};
auto const update_jacobian = [&](LocalJacobianMatrix& jacobian) {
calculateJacobianBurgers(
t, x, dt, jacobian, sig_eff, sigd_j, state.eps_K_j,
- state); // for solution dependent Jacobians
+ local_lubby2_properties); // for solution dependent Jacobians
};
auto const update_solution = [&](LocalResidualVector const& increment) {
@@ -96,7 +101,7 @@ bool Lubby2<DisplacementDim>::computeConstitutiveRelation(
// Calculate effective stress and update material properties
sig_eff = MaterialLib::SolidModels::Invariants<
KelvinVectorSize>::equivalentStress(sigd_j);
- updateBurgersProperties(t, x, sig_eff * state.GM, state);
+ local_lubby2_properties.update(sig_eff);
};
// TODO Make the following choice of maximum iterations and convergence
@@ -127,7 +132,8 @@ bool Lubby2<DisplacementDim>::computeConstitutiveRelation(
double const eps_i_trace = Invariants::trace(eps);
sigma.noalias() =
- state.GM * sigd_j + state.KM * eps_i_trace * Invariants::identity2;
+ local_lubby2_properties.GM0 * sigd_j +
+ local_lubby2_properties.KM0 * eps_i_trace * Invariants::identity2;
// Calculate dGdE for time step
Eigen::Matrix<double, KelvinVectorSize * 3, KelvinVectorSize,
@@ -140,25 +146,12 @@ bool Lubby2<DisplacementDim>::computeConstitutiveRelation(
.template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
auto const& P_sph = Invariants::spherical_projection;
- C.noalias() = state.GM * dzdE * P_dev + 3. * state.KM * P_sph;
+ C.noalias() = local_lubby2_properties.GM0 * dzdE * P_dev +
+ 3. * local_lubby2_properties.KM0 * P_sph;
return true;
}
-template <int DisplacementDim>
-void Lubby2<DisplacementDim>::updateBurgersProperties(
- double const t,
- ProcessLib::SpatialPosition const& x,
- double s_eff,
- MaterialStateVariables& state)
-{
- state.GM = _mp.GM0(t, x)[0];
- state.KM = _mp.KM0(t, x)[0];
- state.GK = _mp.GK0(t, x)[0] * std::exp(_mp.mK(t, x)[0] * s_eff);
- state.etaK = _mp.etaK0(t, x)[0] * std::exp(_mp.mvK(t, x)[0] * s_eff);
- state.etaM = _mp.etaM0(t, x)[0] * std::exp(_mp.mvM(t, x)[0] * s_eff);
-}
-
template <int DisplacementDim>
void Lubby2<DisplacementDim>::calculateResidualBurgers(
const double dt,
@@ -169,7 +162,7 @@ void Lubby2<DisplacementDim>::calculateResidualBurgers(
KelvinVector& strain_Max_curr,
const KelvinVector& strain_Max_t,
ResidualVector& res,
- MaterialStateVariables const& state)
+ detail::LocalLubby2Properties<DisplacementDim> const& properties)
{
// calculate stress residual
res.template block<KelvinVectorSize, 1>(0, 0).noalias() =
@@ -178,13 +171,13 @@ void Lubby2<DisplacementDim>::calculateResidualBurgers(
// calculate Kelvin strain residual
res.template block<KelvinVectorSize, 1>(KelvinVectorSize, 0).noalias() =
1. / dt * (strain_Kel_curr - strain_Kel_t) -
- 1. / (2. * state.etaK) *
- (state.GM * stress_curr - 2. * state.GK * strain_Kel_curr);
+ 1. / (2. * properties.etaK) * (properties.GM0 * stress_curr -
+ 2. * properties.GK * strain_Kel_curr);
// calculate Maxwell strain residual
res.template block<KelvinVectorSize, 1>(2 * KelvinVectorSize, 0).noalias() =
1. / dt * (strain_Max_curr - strain_Max_t) -
- 0.5 * state.GM / state.etaM * stress_curr;
+ 0.5 * properties.GM0 / properties.etaM * stress_curr;
}
template <int DisplacementDim>
@@ -196,7 +189,7 @@ void Lubby2<DisplacementDim>::calculateJacobianBurgers(
double s_eff,
const KelvinVector& sig_i,
const KelvinVector& eps_K_i,
- MaterialStateVariables const& state)
+ detail::LocalLubby2Properties<DisplacementDim> const& properties)
{
Jac.setZero();
@@ -216,42 +209,45 @@ void Lubby2<DisplacementDim>::calculateJacobianBurgers(
// build G_21
Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
- .noalias() = -0.5 * state.GM / state.etaK * KelvinMatrix::Identity();
+ .noalias() =
+ -0.5 * properties.GM0 / properties.etaK * KelvinMatrix::Identity();
if (s_eff > 0.)
{
KelvinVector const eps_K_aid =
- 1. / (state.etaK * state.etaK) *
- (state.GM * sig_i - 2. * state.GK * eps_K_i);
+ 1. / (properties.etaK * properties.etaK) *
+ (properties.GM0 * sig_i - 2. * properties.GK * eps_K_i);
- KelvinVector const dG_K =
- 1.5 * _mp.mK(t, x)[0] * state.GK * state.GM / s_eff * sig_i;
- KelvinVector const dmu_vK =
- 1.5 * _mp.mvK(t, x)[0] * state.GM * state.etaK / s_eff * sig_i;
+ KelvinVector const dG_K = 1.5 * _mp.mK(t, x)[0] * properties.GK *
+ properties.GM0 / s_eff * sig_i;
+ KelvinVector const dmu_vK = 1.5 * _mp.mvK(t, x)[0] * properties.GM0 *
+ properties.etaK / s_eff * sig_i;
Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
0)
.noalias() += 0.5 * eps_K_aid * dmu_vK.transpose() +
- 1. / state.etaK * eps_K_i * dG_K.transpose();
+ 1. / properties.etaK * eps_K_i * dG_K.transpose();
}
// build G_22
Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
KelvinVectorSize)
.diagonal()
- .setConstant(1. / dt + state.GK / state.etaK);
+ .setConstant(1. / dt + properties.GK / properties.etaK);
// nothing to do for G_23
// build G_31
Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
0)
- .noalias() = -0.5 * state.GM / state.etaM * KelvinMatrix::Identity();
+ .noalias() =
+ -0.5 * properties.GM0 / properties.etaM * KelvinMatrix::Identity();
if (s_eff > 0.)
{
- KelvinVector const dmu_vM =
- 1.5 * _mp.mvM(t, x)[0] * state.GM * state.etaM / s_eff * sig_i;
+ KelvinVector const dmu_vM = 1.5 * _mp.mvM(t, x)[0] * properties.GM0 *
+ properties.etaM / s_eff * sig_i;
Jac.template block<KelvinVectorSize, KelvinVectorSize>(
2 * KelvinVectorSize, 0)
- .noalias() += 0.5 * state.GM / (state.etaM * state.etaM) * sig_i *
+ .noalias() += 0.5 * properties.GM0 /
+ (properties.etaM * properties.etaM) * sig_i *
dmu_vM.transpose();
}
diff --git a/MaterialLib/SolidModels/Lubby2.h b/MaterialLib/SolidModels/Lubby2.h
index f216d5f9f7335d8c52855837fd824c2b4f740351..cc39a6995a432c52c9d3cc7b2d31e048e5ec1c38 100644
--- a/MaterialLib/SolidModels/Lubby2.h
+++ b/MaterialLib/SolidModels/Lubby2.h
@@ -56,46 +56,52 @@ struct Lubby2MaterialProperties
P const& mvM;
};
+namespace detail
+{
template <int DisplacementDim>
-class Lubby2 final : public MechanicsBase<DisplacementDim>
+struct LocalLubby2Properties
{
-public:
- //
- // Variables specific to the material model.
- //
- struct MaterialProperties
+ LocalLubby2Properties(double const t,
+ ProcessLib::SpatialPosition const& x,
+ Lubby2MaterialProperties const& mp)
+ : GM0(mp.GM0(t, x)[0]),
+ KM0(mp.KM0(t, x)[0]),
+ GK0(mp.GK0(t, x)[0]),
+ etaK0(mp.etaK0(t, x)[0]),
+ etaM0(mp.etaM0(t, x)[0]),
+ mK(mp.mK(t, x)[0]),
+ mvK(mp.mvK(t, x)[0]),
+ mvM(mp.mvM(t, x)[0])
{
- using P = ProcessLib::Parameter<double>;
- MaterialProperties(P const& GK0_,
- P const& GM0_,
- P const& KM0_,
- P const& etaK0_,
- P const& etaM0_,
- P const& mK_,
- P const& mvK_,
- P const& mvM_)
- : GK0(GK0_),
- GM0(GM0_),
- KM0(KM0_),
- etaK0(etaK0_),
- etaM0(etaM0_),
- mK(mK_),
- mvK(mvK_),
- mvM(mvM_)
- {
- }
+ }
- // basic material parameters
- P const& GK0;
- P const& GM0;
- P const& KM0;
- P const& etaK0;
- P const& etaM0;
- P const& mK;
- P const& mvK;
- P const& mvM;
- };
+ void update(double const s_eff)
+ {
+ GK = GK0 * std::exp(mK * GM0 * s_eff);
+ etaK = etaK0 * std::exp(mvK * GM0 * s_eff);
+ etaM = etaM0 * std::exp(mvM * GM0 * s_eff);
+ }
+ double const GM0;
+ double const KM0;
+ double const GK0;
+ double const etaK0;
+ double const etaM0;
+ double const mK;
+ double const mvK;
+ double const mvM;
+
+ // Solution dependent values.
+ double GK;
+ double etaK;
+ double etaM;
+};
+}
+
+template <int DisplacementDim>
+class Lubby2 final : public MechanicsBase<DisplacementDim>
+{
+public:
struct MaterialStateVariables
: public MechanicsBase<DisplacementDim>::MaterialStateVariables
{
@@ -107,6 +113,12 @@ public:
eps_M_j.resize(KelvinVectorSize);
}
+ void loadState()
+ {
+ eps_K_j = eps_K_t;
+ eps_M_j = eps_M_t;
+ }
+
void pushBackState() override
{
eps_K_t = eps_K_j;
@@ -123,13 +135,6 @@ public:
/// iteration
KelvinVector eps_M_t;
KelvinVector eps_M_j;
-
- // solution dependent values
- double GM;
- double KM;
- double GK;
- double etaK;
- double etaM;
};
std::unique_ptr<
@@ -175,32 +180,28 @@ public:
material_state_variables) override;
private:
- /// Updates Burgers material parameters in LUBBY2 fashion.
- void updateBurgersProperties(double const t,
- ProcessLib::SpatialPosition const& x,
- double const s_eff,
- MaterialStateVariables& _state);
-
/// Calculates the 18x1 residual vector.
- void calculateResidualBurgers(double const dt,
- const KelvinVector& strain_curr,
- const KelvinVector& stress_curr,
- KelvinVector& strain_Kel_curr,
- const KelvinVector& strain_Kel_t,
- KelvinVector& strain_Max_curr,
- const KelvinVector& strain_Max_t,
- ResidualVector& res,
- MaterialStateVariables const& _state);
+ void calculateResidualBurgers(
+ double const dt,
+ const KelvinVector& strain_curr,
+ const KelvinVector& stress_curr,
+ KelvinVector& strain_Kel_curr,
+ const KelvinVector& strain_Kel_t,
+ KelvinVector& strain_Max_curr,
+ const KelvinVector& strain_Max_t,
+ ResidualVector& res,
+ detail::LocalLubby2Properties<DisplacementDim> const& properties);
/// Calculates the 18x18 Jacobian.
- void calculateJacobianBurgers(double const t,
- ProcessLib::SpatialPosition const& x,
- double const dt,
- JacobianMatrix& Jac,
- double s_eff,
- const KelvinVector& sig_i,
- const KelvinVector& eps_K_i,
- MaterialStateVariables const& _state);
+ void calculateJacobianBurgers(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const dt,
+ JacobianMatrix& Jac,
+ double s_eff,
+ const KelvinVector& sig_i,
+ const KelvinVector& eps_K_i,
+ detail::LocalLubby2Properties<DisplacementDim> const& properties);
/// Calculates the 18x6 derivative of the residuals with respect to total
/// strain.