Merge pull request #1556 from endJunction/Ehlers

Ehlers Single surface yield function material model.

Applications/CLI/Tests.cmake

@@ -486,6 +486,47 @@ if(NOT OGS_USE_MPI)
         cube_1e3_expected_pcs_0_ts_101_t_1.000000.vtu cube_1e3_pcs_0_ts_101_t_1.000000.vtu displacement displacement
     )
 
+    # Mechanics; Small deformations, Ehlers (SDE)
+    AddTest(
+        NAME Mechanics_SDE_cube_1e0
+        PATH Mechanics/Ehlers
+        EXECUTABLE ogs
+        EXECUTABLE_ARGS cube_1e0.prj
+        TESTER vtkdiff
+        ABSTOL 1e-15 RELTOL 1e-13
+        DIFF_DATA
+        expected_cube_1e0_pcs_0_ts_101_t_2.550000.vtu cube_1e0_pcs_0_ts_101_t_2.550000.vtu displacement displacement
+        expected_cube_1e0_pcs_0_ts_101_t_2.550000.vtu cube_1e0_pcs_0_ts_101_t_2.550000.vtu sigma_xx sigma_xx
+        expected_cube_1e0_pcs_0_ts_203_t_5.100000.vtu cube_1e0_pcs_0_ts_203_t_5.100000.vtu displacement displacement
+        expected_cube_1e0_pcs_0_ts_203_t_5.100000.vtu cube_1e0_pcs_0_ts_203_t_5.100000.vtu sigma_xx sigma_xx
+    )
+    AddTest(
+        NAME Mechanics_SDE_cube_1e1
+        PATH Mechanics/Ehlers
+        EXECUTABLE ogs
+        EXECUTABLE_ARGS cube_1e1.prj
+        TESTER vtkdiff
+        ABSTOL 1e-15 RELTOL 1e-13
+        DIFF_DATA
+        expected_cube_1e1_pcs_0_ts_101_t_2.550000.vtu cube_1e1_pcs_0_ts_101_t_2.550000.vtu displacement displacement
+        expected_cube_1e1_pcs_0_ts_101_t_2.550000.vtu cube_1e1_pcs_0_ts_101_t_2.550000.vtu sigma_xx sigma_xx
+        expected_cube_1e1_pcs_0_ts_203_t_5.100000.vtu cube_1e1_pcs_0_ts_203_t_5.100000.vtu displacement displacement
+        expected_cube_1e1_pcs_0_ts_203_t_5.100000.vtu cube_1e1_pcs_0_ts_203_t_5.100000.vtu sigma_xx sigma_xx
+    )
+    AddTest(
+        NAME LARGE_Mechanics_SDE_cube_1e3
+        PATH Mechanics/Ehlers
+        EXECUTABLE ogs
+        EXECUTABLE_ARGS cube_1e3.prj
+        TESTER vtkdiff
+        ABSTOL 1e-15 RELTOL 1e-13
+        DIFF_DATA
+        expected_cube_1e3_pcs_0_ts_101_t_2.550000.vtu cube_1e3_pcs_0_ts_101_t_2.550000.vtu displacement displacement
+        expected_cube_1e3_pcs_0_ts_101_t_2.550000.vtu cube_1e3_pcs_0_ts_101_t_2.550000.vtu sigma_xx sigma_xx
+        expected_cube_1e3_pcs_0_ts_203_t_5.100000.vtu cube_1e3_pcs_0_ts_203_t_5.100000.vtu displacement displacement
+        expected_cube_1e3_pcs_0_ts_203_t_5.100000.vtu cube_1e3_pcs_0_ts_203_t_5.100000.vtu sigma_xx sigma_xx
+    )
+
     # SQUARE 1x1 GROUNDWATER FLOW TEST -- AXIALLY SYMMETRIC
     # test results are compared to 3D simulation on a wedge-shaped domain
     AddTest(

Documentation/bibliography.bib

@@ -2,6 +2,17 @@
 % Encoding: UTF-8
 
 
+@Article{Ehlers1995,
+  Title={A single-surface yield function for geomaterials},
+  Author={Ehlers, W},
+  Journal={Archive of Applied Mechanics},
+  Volume={65},
+  Number={4},
+  Pages={246--259},
+  Year={1995},
+  Publisher={Springer}
+}
+
 @Book{Ericson:2004:RCD:1121584,
   Title                    = {Real-Time Collision Detection (The Morgan Kaufmann Series in Interactive 3-D Technology) (The Morgan Kaufmann Series in Interactive 3D Technology)},
   Author                   = {Ericson, Christer},

MaterialLib/SolidModels/CreateEhlers.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_CREATEEHLERS_H_
+#define MATERIALLIB_SOLIDMODELS_CREATEEHLERS_H_
+
+#include <logog/include/logog.hpp>
+
+#include "MechanicsBase.h"
+#include "ProcessLib/Utils/ProcessUtils.h"  // required for findParameter
+#include "Ehlers.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+std::unique_ptr<MechanicsBase<DisplacementDim>> createEhlers(
+    std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+    BaseLib::ConfigTree const& config)
+{
+    config.checkConfigParameter("type", "Ehlers");
+    DBUG("Create Ehlers material");
+
+    auto& shear_modulus = ProcessLib::findParameter<double>(
+        config, "shear_modulus", parameters, 1);
+
+    DBUG("Use \'%s\' as shear modulus parameter.", shear_modulus.name.c_str());
+
+    auto& bulk_modulus = ProcessLib::findParameter<double>(
+        config, "bulk_modulus", parameters, 1);
+
+    DBUG("Use \'%s\' as bulk modulus parameter.", bulk_modulus.name.c_str());
+
+    auto& kappa = ProcessLib::findParameter<double>(
+        config, "kappa", parameters, 1);
+
+    DBUG("Use \'%s\' as kappa.",
+         kappa.name.c_str());
+
+    auto& beta = ProcessLib::findParameter<double>(
+        config, "beta", parameters, 1);
+
+    DBUG("Use \'%s\' as beta.",
+         beta.name.c_str());
+
+    auto& gamma = ProcessLib::findParameter<double>(
+        config, "gamma", parameters, 1);
+
+    DBUG("Use \'%s\' as gamma.",
+         gamma.name.c_str());
+
+    auto& hardening_modulus = ProcessLib::findParameter<double>(
+        config, "hardening_modulus", parameters, 1);
+
+    DBUG("Use \'%s\' as hardening modulus parameter.",
+         hardening_modulus.name.c_str());
+
+    auto& alpha = ProcessLib::findParameter<double>(
+        config, "alpha", parameters, 1);
+
+    DBUG("Use \'%s\' as alpha.",
+         alpha.name.c_str());
+
+    auto& delta = ProcessLib::findParameter<double>(
+        config, "delta", parameters, 1);
+
+    DBUG("Use \'%s\' as delta.",
+         delta.name.c_str());
+
+    auto& eps = ProcessLib::findParameter<double>(
+        config, "eps", parameters, 1);
+
+    DBUG("Use \'%s\' as eps.",
+         eps.name.c_str());
+
+    auto& m = ProcessLib::findParameter<double>(
+        config, "m", parameters, 1);
+
+    DBUG("Use \'%s\' as m.",
+         m.name.c_str());
+
+    auto& alphap = ProcessLib::findParameter<double>(
+        config, "alphap", parameters, 1);
+
+    DBUG("Use \'%s\' as alphap.",
+         alphap.name.c_str());
+
+    auto& deltap = ProcessLib::findParameter<double>(
+        config, "deltap", parameters, 1);
+
+    DBUG("Use \'%s\' as deltap.",
+         deltap.name.c_str());
+
+    auto& epsp = ProcessLib::findParameter<double>(
+        config, "epsp", parameters, 1);
+
+    DBUG("Use \'%s\' as epsp.",
+         epsp.name.c_str());
+
+    auto& paremeter_mp = ProcessLib::findParameter<double>(
+        config, "mp", parameters, 1);
+
+    DBUG("Use \'%s\' as mp.",
+         paremeter_mp.name.c_str());
+
+    auto& betap = ProcessLib::findParameter<double>(
+        config, "betap", parameters, 1);
+
+    DBUG("Use \'%s\' as betap.",
+         betap.name.c_str());
+
+    auto& gammap = ProcessLib::findParameter<double>(
+        config, "gammap", parameters, 1);
+
+    DBUG("Use \'%s\' as gammap.",
+         gammap.name.c_str());
+
+    typename SolidEhlers<DisplacementDim>::MaterialProperties mp{
+        shear_modulus, bulk_modulus, alpha,  beta,
+        gamma,         delta,        eps,    m,
+        alphap,        betap,        gammap, deltap,
+        epsp,          paremeter_mp, kappa,  hardening_modulus};
+
+    return std::unique_ptr<MechanicsBase<DisplacementDim>>{
+        new SolidEhlers<DisplacementDim>{mp}};
+}
+
+}  // namespace Solids
+}  // namespace MaterialLib
+
+#endif  // MATERIALLIB_SOLIDMODELS_CREATEEHLERS_H_

MaterialLib/SolidModels/Ehlers.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 "Ehlers.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template class SolidEhlers<2>;
+template class SolidEhlers<3>;
+
+}  // namespace Solids
+}  // namespace MaterialLib

MaterialLib/SolidModels/Ehlers.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
+ *
+ */
+
+/*
+ * Implementation of Ehler's single-surface model.
+ * see Ehler's paper "A single-surface yield function for geomaterials" for more
+ * details. \cite{Ehlers1995}
+ *
+ * Refer to "Single-surface benchmark of OpenGeoSys documentation
+ * (https://docs.opengeosys.org/docs/benchmarks/small-deformations/mechanics-plasticity-single-surface)"
+ * for more details for the tests.
+ */
+
+#ifndef MATERIALLIB_SOLIDMODELS_EHLERS_H_
+#define MATERIALLIB_SOLIDMODELS_EHLERS_H_
+
+#include <cfloat>
+#include <memory>
+#ifndef NDEBUG
+#include <ostream>
+#endif
+
+#include <Eigen/Dense>
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/Error.h"
+#include "KelvinVector.h"
+#include "MechanicsBase.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+
+template <int DisplacementDim>
+class SolidEhlers final : public MechanicsBase<DisplacementDim>
+{
+public:
+    static int const KelvinVectorSize =
+        ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
+    static int const JacobianResidualSize =
+        2 * KelvinVectorSize + 3;  // 2 is the number of components in the
+                                   // jacobian/residual, not the space
+                                   // dimension. And 3 is for additional
+                                   // variables.
+    using ResidualVectorType = Eigen::Matrix<double, JacobianResidualSize, 1>;
+    using JacobianMatrix = Eigen::Matrix<double, JacobianResidualSize,
+                                         JacobianResidualSize, Eigen::RowMajor>;
+
+public:
+    //
+    // Variables specific to the material model.
+    //
+    struct MaterialProperties
+    {
+        using P = ProcessLib::Parameter<double>;
+
+        using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+        using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
+
+        /// material parameters in relation to Ehler's single-surface model
+        /// see Ehler's paper "A single-surface yield function for geomaterials"
+        /// for more details.
+        MaterialProperties(P const& G_, P const& K_, P const& alpha_,
+                           P const& beta_, P const& gamma_, P const& delta_,
+                           P const& epsilon_, P const& m_, P const& alpha_p_,
+                           P const& beta_p_, P const& gamma_p_,
+                           P const& delta_p_, P const& epsilon_p_,
+                           P const& m_p_, P const& kappa_,
+                           P const& hardening_coefficient_)
+            : G(G_),
+              K(K_),
+              alpha(alpha_),
+              beta(beta_),
+              gamma(gamma_),
+              delta(delta_),
+              epsilon(epsilon_),
+              m(m_),
+              alpha_p(alpha_p_),
+              beta_p(beta_p_),
+              gamma_p(gamma_p_),
+              delta_p(delta_p_),
+              epsilon_p(epsilon_p_),
+              m_p(m_p_),
+              kappa(kappa_),
+              hardening_coefficient(hardening_coefficient_)
+        {
+        }
+        // basic material parameters
+        P const& G;  ///< shear modulus
+        P const& K;  ///< bulk modulus
+
+        P const& alpha;
+        P const& beta;
+        P const& gamma;
+        P const& delta;
+        P const& epsilon;
+        P const& m;
+
+        P const& alpha_p;
+        P const& beta_p;
+        P const& gamma_p;
+        P const& delta_p;
+        P const& epsilon_p;
+        P const& m_p;
+
+        P const& kappa;
+        P const& hardening_coefficient;
+        // Drucker-Prager: Import kappa and beta in terms of Drucker-Prager
+        // criterion solution dependent values
+        double k;
+
+        void calculateIsotropicHardening(double const t,
+                                         ProcessLib::SpatialPosition const& x,
+                                         const double e_pv_curr);
+    };
+
+    struct MaterialStateVariables
+        : public MechanicsBase<DisplacementDim>::MaterialStateVariables
+    {
+        MaterialStateVariables()
+            : eps_p_D(KelvinVector::Zero()), eps_p_D_prev(KelvinVector::Zero())
+        {
+        }
+
+        void loadState()
+        {
+            eps_p_D = eps_p_D_prev;
+            eps_p_V = eps_p_V_prev;
+            eps_p_eff = eps_p_eff_prev;
+            lambda = 0;
+        }
+
+        void pushBackState() override
+        {
+            eps_p_D_prev = eps_p_D;
+            eps_p_V_prev = eps_p_V;
+            eps_p_eff_prev = eps_p_eff;  // effective part of trace(eps_p)
+            lambda = 0;
+        }
+
+        using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+
+        KelvinVector eps_p_D;  ///< deviatoric plastic strain
+        double eps_p_V = 0;    ///< volumetric strain
+        double eps_p_eff = 0;  ///< effective plastic strain
+
+        // Initial values from previous timestep
+        KelvinVector eps_p_D_prev;  ///< \copydoc eps_p_D
+        double eps_p_V_prev = 0;    ///< \copydoc eps_p_V
+        double eps_p_eff_prev = 0;  ///< \copydoc eps_p_eff
+
+        double lambda = 0;  ///< plastic multiplier
+
+#ifndef NDEBUG
+        friend std::ostream& operator<<(std::ostream& os,
+                                        MaterialStateVariables const& m)
+        {
+            os << "State:\n"
+               << "eps_p_D: " << m.eps_p_D << "\n"
+               << "e_p: " << m.e_p << "\n"
+               << "eps_p_eff: " << m.eps_p_eff << "\n"
+               << "eps_p_D_prev: " << m.eps_p_D_prev << "\n"
+               << "e_p_prev: " << m.e_p_prev << "\n"
+               << "eps_p_eff_prev: " << m.eps_p_eff_prev << "\n"
+               << "lambda: " << m.lambda << "\n";
+            return os;
+        }
+#endif  // NDEBUG
+    };
+
+    std::unique_ptr<
+        typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
+    createMaterialStateVariables() override
+    {
+        return std::unique_ptr<
+            typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
+            new MaterialStateVariables};
+    }
+
+public:
+    using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+    using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
+
+public:
+    explicit SolidEhlers(MaterialProperties const& material_properties)
+        : _mp(material_properties)
+    {
+    }
+
+    bool 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) override;
+
+private:
+    MaterialProperties _mp;
+};
+
+}  // namespace Solids
+}  // namespace MaterialLib
+#include "Ehlers-impl.h"
+
+#endif  // MATERIALLIB_SOLIDMODELS_EHLERS_H_

ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp

@@ -13,6 +13,7 @@
 
 #include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
 #include "MaterialLib/SolidModels/CreateLubby2.h"
+#include "MaterialLib/SolidModels/CreateEhlers.h"
 #include "ProcessLib/Utils/ParseSecondaryVariables.h"
 
 #include "SmallDeformationProcess.h"
@@ -74,7 +75,12 @@ createSmallDeformationProcess(
 
     std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
         material = nullptr;
-    if (type == "LinearElasticIsotropic")
+    if (type == "Ehlers")
+    {
+        material = MaterialLib::Solids::createEhlers<DisplacementDim>(
+            parameters, constitutive_relation_config);
+    }
+    else if (type == "LinearElasticIsotropic")
     {
         material =
             MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(

ProcessLib/SmallDeformation/SmallDeformationFEM.h

@@ -112,6 +112,24 @@ struct SmallDeformationLocalAssemblerInterface
 
     virtual std::vector<double> const& getIntPtSigmaYZ(
         std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonXX(
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonYY(
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonZZ(
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonXY(
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonXZ(
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilonYZ(
+        std::vector<double>& cache) const = 0;
 };
 
 template <typename ShapeFunction, typename IntegrationMethod,
@@ -318,6 +336,44 @@ public:
         return getIntPtSigma(cache, 5);
     }
 
+    std::vector<double> const& getIntPtEpsilonXX(
+        std::vector<double>& cache) const override
+    {
+        return getIntPtEpsilon(cache, 0);
+    }
+
+    std::vector<double> const& getIntPtEpsilonYY(
+        std::vector<double>& cache) const override
+    {
+        return getIntPtEpsilon(cache, 1);
+    }
+
+    std::vector<double> const& getIntPtEpsilonZZ(
+        std::vector<double>& cache) const override
+    {
+        return getIntPtEpsilon(cache, 2);
+    }
+
+    std::vector<double> const& getIntPtEpsilonXY(
+        std::vector<double>& cache) const override
+    {
+        return getIntPtEpsilon(cache, 3);
+    }
+
+    std::vector<double> const& getIntPtEpsilonXZ(
+        std::vector<double>& cache) const override
+    {
+        assert(DisplacementDim == 3);
+        return getIntPtEpsilon(cache, 4);
+    }
+
+    std::vector<double> const& getIntPtEpsilonYZ(
+        std::vector<double>& cache) const override
+    {
+        assert(DisplacementDim == 3);
+        return getIntPtEpsilon(cache, 5);
+    }
+
 private:
     std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
                                              std::size_t const component) const
@@ -335,6 +391,20 @@ private:
         return cache;
     }
 
+    std::vector<double> const& getIntPtEpsilon(
+        std::vector<double>& cache, std::size_t const component) const
+    {
+        cache.clear();
+        cache.reserve(_ip_data.size());
+
+        for (auto const& ip_data : _ip_data) {
+            cache.push_back(ip_data._eps[component]);
+        }
+
+        return cache;
+    }
+
+
     SmallDeformationProcessData<DisplacementDim>& _process_data;
 
     std::vector<IntegrationPointData<BMatricesType, DisplacementDim>> _ip_data;

ProcessLib/SmallDeformation/SmallDeformationProcess.h

@@ -116,6 +116,30 @@ private:
                     getExtrapolator(), _local_assemblers,
                     &SmallDeformationLocalAssemblerInterface::getIntPtSigmaYZ));
         }
+
+        Base::_secondary_variables.addSecondaryVariable(
+            "epsilon_xx", 1,
+            makeExtrapolator(
+                getExtrapolator(), _local_assemblers,
+                &SmallDeformationLocalAssemblerInterface::getIntPtEpsilonXX));
+
+        Base::_secondary_variables.addSecondaryVariable(
+            "epsilon_yy", 1,
+            makeExtrapolator(
+                getExtrapolator(), _local_assemblers,
+                &SmallDeformationLocalAssemblerInterface::getIntPtEpsilonYY));
+
+        Base::_secondary_variables.addSecondaryVariable(
+            "epsilon_zz", 1,
+            makeExtrapolator(
+                getExtrapolator(), _local_assemblers,
+                &SmallDeformationLocalAssemblerInterface::getIntPtEpsilonZZ));
+
+        Base::_secondary_variables.addSecondaryVariable(
+            "epsilon_xy", 1,
+            makeExtrapolator(
+                getExtrapolator(), _local_assemblers,
+                &SmallDeformationLocalAssemblerInterface::getIntPtEpsilonXY));
     }
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,

Data @ 74db8c48

-Subproject commit 6ffd7ec67ebd30d3c83754cf3127fd885e255528
+Subproject commit 74db8c483fcdd3231b40215a78a9357865a80f5e