diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/c_Lubby2.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/c_Lubby2.md
new file mode 100644
index 0000000000000000000000000000000000000000..5fe806062936d160f3c43f3f7ceaf18c708b644e
--- /dev/null
+++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/c_Lubby2.md
@@ -0,0 +1,3 @@
+Visco elastic material model based on Burgers' rheological model. Viscosities
+and Kelvin shear modulus are stress dependent according to LUBBY2
+model. Heusermann 1983.
diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mk.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mk.md
new file mode 100644
index 0000000000000000000000000000000000000000..08aa5f31811545b9464b045b1d71ee2d2914f9b1
--- /dev/null
+++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mk.md
@@ -0,0 +1 @@
+Parameter characterizing the stress dependence of the Kelvin shear modulus.
diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mvk.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mvk.md
new file mode 100644
index 0000000000000000000000000000000000000000..7890d68f729fb41a0bf7f1e17b51e759ee6edab0
--- /dev/null
+++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mvk.md
@@ -0,0 +1 @@
+Parameter characterizing the stress dependence of the Kelvin shear viscosity.
diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mvm.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mvm.md
new file mode 100644
index 0000000000000000000000000000000000000000..a2efbaee8d68602eae848415a6ad3d85b36559d8
--- /dev/null
+++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mvm.md
@@ -0,0 +1 @@
+Parameter characterizing the stress dependence of the Maxwell shear viscosity.
diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_kelvin_shear_modulus.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_kelvin_shear_modulus.md
new file mode 100644
index 0000000000000000000000000000000000000000..48367a0f36fa5db059fadc30c6f3dd92e11db3e7
--- /dev/null
+++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_kelvin_shear_modulus.md
@@ -0,0 +1 @@
+Initial shear modulus of the Kelvin element in the Lubby2 model.
diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_kelvin_viscosity.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_kelvin_viscosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..e8577b7d43d73b1f4169301c01a1ca0cc2cc619d
--- /dev/null
+++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_kelvin_viscosity.md
@@ -0,0 +1 @@
+Initial shear viscosity of the Kelvin element in the Lubby2 model.
diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_bulk_modulus.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_bulk_modulus.md
new file mode 100644
index 0000000000000000000000000000000000000000..722f42409620670a05b9a4471b7fffc5acd76190
--- /dev/null
+++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_bulk_modulus.md
@@ -0,0 +1 @@
+Bulk modulus of the Maxwell element in the Lubby2 model.
diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_shear_modulus.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_shear_modulus.md
new file mode 100644
index 0000000000000000000000000000000000000000..e67fcd4b0c30c3b4c1d2862201dd89fbac1c7665
--- /dev/null
+++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_shear_modulus.md
@@ -0,0 +1 @@
+Shear modulus of the Maxwell element in the Lubby2 model.
diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_viscosity.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_viscosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..6f1701b5cf6a76a759e152ffce9c044e003dbfe0
--- /dev/null
+++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_viscosity.md
@@ -0,0 +1 @@
+Initial shear viscosity of the Maxwell element in the Lubby2 model.
diff --git a/MaterialLib/SolidModels/CreateLubby2.h b/MaterialLib/SolidModels/CreateLubby2.h
new file mode 100644
index 0000000000000000000000000000000000000000..d06dad3c2e3be9fa054c65cc79332fa0afd516c4
--- /dev/null
+++ b/MaterialLib/SolidModels/CreateLubby2.h
@@ -0,0 +1,101 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALLIB_SOLIDMODELS_CREATELUBBY2_H_
+#define MATERIALLIB_SOLIDMODELS_CREATELUBBY2_H_
+
+#include <logog/include/logog.hpp>
+
+#include "Lubby2.h"
+#include "MechanicsBase.h"
+#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+std::unique_ptr<MechanicsBase<DisplacementDim>> createLubby2(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config)
+{
+ config.checkConfigParameter("type", "Lubby2");
+ DBUG("Create Lubby2 material");
+
+ // Kelvin shear modulus.
+ auto& kelvin_shear_modulus = ProcessLib::findParameter<double>(
+ config, "kelvin_shear_modulus", parameters, 1);
+
+ DBUG("Use '%s' as kelvin shear modulus parameter.",
+ kelvin_shear_modulus.name.c_str());
+
+ // Kelvin viscosity.
+ auto& kelvin_viscosity = ProcessLib::findParameter<double>(
+ config, "kelvin_viscosity", parameters, 1);
+
+ DBUG("Use '%s' as kelvin viscosity parameter.",
+ kelvin_viscosity.name.c_str());
+
+ // Maxwell shear modulus.
+ auto& maxwell_shear_modulus = ProcessLib::findParameter<double>(
+ config, "maxwell_shear_modulus", parameters, 1);
+
+ DBUG("Use '%s' as maxwell shear modulus parameter.",
+ maxwell_shear_modulus.name.c_str());
+
+ // Maxwell bulk modulus.
+ auto& maxwell_bulk_modulus = ProcessLib::findParameter<double>(
+ config, "maxwell_bulk_modulus", parameters, 1);
+
+ DBUG("Use '%s' as maxwell bulk modulus parameter.",
+ maxwell_bulk_modulus.name.c_str());
+
+ // Maxwell viscosity.
+ auto& maxwell_viscosity = ProcessLib::findParameter<double>(
+ config, "maxwell_viscosity", parameters, 1);
+
+ DBUG("Use '%s' as maxwell viscosity parameter.",
+ maxwell_viscosity.name.c_str());
+
+ // Dependency parameter for mK.
+ auto& dependency_parameter_mK = ProcessLib::findParameter<double>(
+ config, "dependency_parameter_mk", parameters, 1);
+
+ DBUG("Use '%s' as dependency parameter mK.",
+ dependency_parameter_mK.name.c_str());
+
+ // Dependency parameter for mvK.
+ auto& dependency_parameter_mvK = ProcessLib::findParameter<double>(
+ config, "dependency_parameter_mvk", parameters, 1);
+
+ DBUG("Use '%s' as dependency parameter mvK.",
+ dependency_parameter_mvK.name.c_str());
+
+ // Dependency parameter for mvM.
+ auto& dependency_parameter_mvM = ProcessLib::findParameter<double>(
+ config, "dependency_parameter_mvm", parameters, 1);
+
+ DBUG("Use '%s' as dependency parameter mvM.",
+ dependency_parameter_mvM.name.c_str());
+
+
+ typename Lubby2<DisplacementDim>::MaterialProperties mp{
+ kelvin_shear_modulus, maxwell_shear_modulus,
+ maxwell_bulk_modulus, kelvin_viscosity,
+ maxwell_viscosity, dependency_parameter_mK,
+ dependency_parameter_mvK, dependency_parameter_mvM};
+
+ return std::unique_ptr<MechanicsBase<DisplacementDim>>{
+ new Lubby2<DisplacementDim>{mp}};
+}
+
+} // namespace Solids
+} // namespace MaterialLib
+
+#endif // MATERIALLIB_SOLIDMODELS_CREATELUBBY2_H_
diff --git a/MaterialLib/SolidModels/Lubby2-impl.h b/MaterialLib/SolidModels/Lubby2-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..392fd020edec5ff549a75cc6161e98a41335de43
--- /dev/null
+++ b/MaterialLib/SolidModels/Lubby2-impl.h
@@ -0,0 +1,269 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALLIB_SOLIDMODELS_LUBBY2_IMPL_H_
+#define MATERIALLIB_SOLIDMODELS_LUBBY2_IMPL_H_
+
+#include "NewtonRaphson.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+void Lubby2<DisplacementDim>::computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const dt,
+ KelvinVector const& /*eps_prev*/,
+ KelvinVector const& eps,
+ KelvinVector const& /*sigma_prev*/,
+ KelvinVector& sigma,
+ KelvinMatrix& C,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables&
+ material_state_variables)
+{
+ using Invariants = MaterialLib::SolidModels::Invariants<KelvinVectorSize>;
+
+ assert(dynamic_cast<MaterialStateVariables*>(&material_state_variables) !=
+ nullptr);
+ MaterialStateVariables& state =
+ static_cast<MaterialStateVariables&>(material_state_variables);
+ // calculation of deviatoric parts
+ auto const& P_dev = Invariants::deviatoric_projection;
+ KelvinVector const epsd_i = P_dev * eps;
+
+ // initial guess as elastic predictor
+ KelvinVector sigd_j = 2.0 * (epsd_i - state.eps_M_t - state.eps_K_t);
+
+ // Calculate effective stress and update material properties
+ double sig_eff = Invariants::equivalentStress(sigd_j);
+ updateBurgersProperties(t, x, sig_eff * state.GM, state);
+
+ using LocalJacobianMatrix =
+ Eigen::Matrix<double, KelvinVectorSize * 3, KelvinVectorSize * 3,
+ Eigen::RowMajor>;
+
+ // Linear solver for the newton loop is required after the loop with the
+ // same matrix. This saves one decomposition.
+ Eigen::PartialPivLU<LocalJacobianMatrix> linear_solver(KelvinVectorSize *
+ 3);
+
+ // Different solvers are available for the solution of the local system.
+ // TODO Make the following choice of linear solvers available from the
+ // input file configuration:
+ // K_loc.partialPivLu().solve(-res_loc);
+ // K_loc.fullPivLu().solve(-res_loc);
+ // K_loc.householderQr().solve(-res_loc);
+ // K_loc.colPivHouseholderQr().solve(res_loc);
+ // K_loc.fullPivHouseholderQr().solve(-res_loc);
+ // K_loc.llt().solve(-res_loc);
+ // K_loc.ldlt().solve(-res_loc);
+
+ { // Local Newton solver
+ using LocalResidualVector =
+ Eigen::Matrix<double, KelvinVectorSize * 3, 1>;
+
+ 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);
+ };
+
+ auto const update_jacobian = [&](LocalJacobianMatrix& jacobian) {
+ calculateJacobianBurgers(
+ t, x, dt, jacobian, sig_eff, sigd_j, state.eps_K_j,
+ state); // for solution dependent Jacobians
+ };
+
+ auto const update_solution = [&](LocalResidualVector const& increment) {
+ // increment solution vectors
+ sigd_j.noalias() += increment.template block<KelvinVectorSize, 1>(
+ KelvinVectorSize * 0, 0);
+ state.eps_K_j.noalias() +=
+ increment.template block<KelvinVectorSize, 1>(
+ KelvinVectorSize * 1, 0);
+ state.eps_M_j.noalias() +=
+ increment.template block<KelvinVectorSize, 1>(
+ KelvinVectorSize * 2, 0);
+
+ // Calculate effective stress and update material properties
+ sig_eff = MaterialLib::SolidModels::Invariants<
+ KelvinVectorSize>::equivalentStress(sigd_j);
+ updateBurgersProperties(t, x, sig_eff * state.GM, state);
+ };
+
+ // TODO Make the following choice of maximum iterations and convergence
+ // criteria available from the input file configuration:
+ const int maximum_iterations(20);
+ const double tolerance(1.e-10);
+
+ auto newton_solver =
+ NewtonRaphson<decltype(linear_solver), LocalJacobianMatrix,
+ decltype(update_jacobian), LocalResidualVector,
+ decltype(update_residual), decltype(update_solution)>(
+ linear_solver, update_jacobian, update_residual,
+ update_solution, maximum_iterations, tolerance);
+
+ auto const success_iterations = newton_solver.solve(K_loc);
+
+ if (!success_iterations)
+ OGS_FATAL("The local Newton method did not succeed.");
+
+ // If the Newton loop didn't run, the linear solver will not be
+ // initialized.
+ // This happens usually for the first iteration of the first timestep.
+ if (*success_iterations == 0)
+ linear_solver.compute(K_loc);
+ }
+
+ // Hydrostatic part for the stress and the tangent.
+ double const eps_i_trace = Invariants::trace(eps);
+
+ sigma.noalias() =
+ state.GM * sigd_j + state.KM * eps_i_trace * Invariants::identity2;
+
+ // Calculate dGdE for time step
+ Eigen::Matrix<double, KelvinVectorSize * 3, KelvinVectorSize,
+ Eigen::RowMajor> const dGdE = calculatedGdEBurgers();
+
+ // Consistent tangent from local Newton iteration of material
+ // functionals.
+ // Only the upper left block is relevant for the global tangent.
+ auto dzdE = linear_solver.solve(-dGdE)
+ .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;
+}
+
+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,
+ 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)
+{
+ // calculate stress residual
+ res.template block<KelvinVectorSize, 1>(0, 0).noalias() =
+ stress_curr - 2. * (strain_curr - strain_Kel_curr - strain_Max_curr);
+
+ // 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);
+
+ // 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;
+}
+
+template <int DisplacementDim>
+void Lubby2<DisplacementDim>::calculateJacobianBurgers(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ const double dt,
+ JacobianMatrix& Jac,
+ double s_eff,
+ const KelvinVector& sig_i,
+ const KelvinVector& eps_K_i,
+ MaterialStateVariables const& state)
+{
+ Jac.setZero();
+
+ // build G_11
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, 0).setIdentity();
+
+ // build G_12
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
+ .diagonal()
+ .setConstant(2);
+
+ // build G_13
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(0,
+ 2 * KelvinVectorSize)
+ .diagonal()
+ .setConstant(2);
+
+ // build G_21
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
+ .noalias() = -0.5 * state.GM / state.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);
+
+ 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;
+ 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();
+ }
+
+ // build G_22
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
+ KelvinVectorSize)
+ .diagonal()
+ .setConstant(1. / dt + state.GK / state.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();
+ if (s_eff > 0.)
+ {
+ KelvinVector const dmu_vM =
+ 1.5 * _mp.mvM(t, x)[0] * state.GM * state.etaM / s_eff * sig_i;
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(
+ 2 * KelvinVectorSize, 0)
+ .noalias() += 0.5 * state.GM / (state.etaM * state.etaM) * sig_i *
+ dmu_vM.transpose();
+ }
+
+ // nothing to do for G_32
+
+ // build G_33
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
+ 2 * KelvinVectorSize)
+ .diagonal()
+ .setConstant(1. / dt);
+}
+
+} // namespace Solids
+} // namespace MaterialLib
+
+#endif // MATERIALLIB_SOLIDMODELS_LUBBY2_IMPL_H_
diff --git a/MaterialLib/SolidModels/Lubby2.cpp b/MaterialLib/SolidModels/Lubby2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6772150832e6320125681e202dc3ab857565862f
--- /dev/null
+++ b/MaterialLib/SolidModels/Lubby2.cpp
@@ -0,0 +1,20 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "Lubby2.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template class Lubby2<2>;
+template class Lubby2<3>;
+
+} // namespace Solids
+} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/Lubby2.h b/MaterialLib/SolidModels/Lubby2.h
new file mode 100644
index 0000000000000000000000000000000000000000..e2d8d50672e27cf8c3506a8c87add55b6e31da19
--- /dev/null
+++ b/MaterialLib/SolidModels/Lubby2.h
@@ -0,0 +1,205 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALLIB_SOLIDMODELS_LUBBY2_H_
+#define MATERIALLIB_SOLIDMODELS_LUBBY2_H_
+
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/Error.h"
+
+#include "KelvinVector.h"
+#include "MechanicsBase.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+class Lubby2 final : public MechanicsBase<DisplacementDim>
+{
+public:
+ //
+ // Variables specific to the material model.
+ //
+ struct MaterialProperties
+ {
+ 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;
+ };
+
+ struct MaterialStateVariables
+ : public MechanicsBase<DisplacementDim>::MaterialStateVariables
+ {
+ MaterialStateVariables()
+ {
+ eps_K_t.resize(KelvinVectorSize);
+ eps_K_j.resize(KelvinVectorSize);
+ eps_M_t.resize(KelvinVectorSize);
+ eps_M_j.resize(KelvinVectorSize);
+ }
+
+ void pushBackState() override
+ {
+ eps_K_t = eps_K_j;
+ eps_M_t = eps_M_j;
+ }
+
+ using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+ using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
+ /// Deviatoric strain in the viscous kelvin element during the current
+ /// iteration
+ KelvinVector eps_K_t;
+ KelvinVector eps_K_j;
+ /// Deviatoric strain in the viscous maxwell element during the current
+ /// 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<
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
+ createMaterialStateVariables() override
+ {
+ return std::unique_ptr<
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
+ new MaterialStateVariables};
+ }
+
+public:
+ static int const KelvinVectorSize =
+ ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
+ using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+ using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
+
+ static int const JacobianResidualSize =
+ 3 * KelvinVectorSize; // Three is the number of components in the
+ // jacobian/residual, not the space dimension.
+ using ResidualVector = Eigen::Matrix<double, JacobianResidualSize, 1>;
+ using JacobianMatrix = Eigen::Matrix<double,
+ JacobianResidualSize,
+ JacobianResidualSize,
+ Eigen::RowMajor>;
+
+public:
+ explicit Lubby2(MaterialProperties& material_properties)
+ : _mp(material_properties)
+ {
+ }
+
+ void computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x_position,
+ double const dt,
+ KelvinVector const& eps_prev,
+ KelvinVector const& eps,
+ KelvinVector const& sigma_prev,
+ KelvinVector& sigma,
+ KelvinMatrix& C,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables&
+ 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);
+
+ /// 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);
+
+ /// Calculates the 18x6 derivative of the residuals with respect to total
+ /// strain.
+ ///
+ /// Function definition can not be moved into implementation because of a
+ /// MSVC compiler errors. See
+ /// http://stackoverflow.com/questions/1484885/strange-vc-compile-error-c2244
+ /// and https://support.microsoft.com/en-us/kb/930198
+ Eigen::
+ Matrix<double, JacobianResidualSize, KelvinVectorSize, Eigen::RowMajor>
+ calculatedGdEBurgers() const
+ {
+ Eigen::Matrix<double,
+ JacobianResidualSize,
+ KelvinVectorSize,
+ Eigen::RowMajor>
+ dGdE;
+ dGdE.setZero();
+ dGdE.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
+ .diagonal()
+ .setConstant(-2);
+ return dGdE;
+ }
+
+private:
+ MaterialProperties _mp;
+};
+
+extern template class Lubby2<2>;
+extern template class Lubby2<3>;
+
+} // namespace Solids
+} // namespace MaterialLib
+
+#include "Lubby2-impl.h"
+
+#endif // MATERIALLIB_SOLIDMODELS_LUBBY2_H_