[Mat] ViscoElasticBurgers model implementation.

Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/c_Lubby2.md

+Visco elastic material model based on Burgers' rheological model. Viscosities
+and Kelvin shear modulus are stress dependent according to LUBBY2
+model. Heusermann 1983.

Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mk.md

+Parameter characterizing the stress dependence of the Kelvin shear modulus.

Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mvk.md

+Parameter characterizing the stress dependence of the Kelvin shear viscosity.

Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_dependency_parameter_mvm.md

+Parameter characterizing the stress dependence of the Maxwell shear viscosity.

Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_kelvin_shear_modulus.md

+Initial shear modulus of the Kelvin element in the Lubby2 model.

Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_kelvin_viscosity.md

+Initial shear viscosity of the Kelvin element in the Lubby2 model.

Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_bulk_modulus.md

+Bulk modulus of the Maxwell element in the Lubby2 model.

Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_shear_modulus.md

+Shear modulus of the Maxwell element in the Lubby2 model.

Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/Lubby2/t_maxwell_viscosity.md

+Initial shear viscosity of the Maxwell element in the Lubby2 model.

MaterialLib/SolidModels/CreateLubby2.h

+/**
+ * \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_

MaterialLib/SolidModels/Lubby2-impl.h

+/**
+ * \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_

MaterialLib/SolidModels/Lubby2.cpp

+/**
+ * \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

MaterialLib/SolidModels/Lubby2.h

+/**
+ * \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_