diff --git a/MaterialLib/CMakeLists.txt b/MaterialLib/CMakeLists.txt
index 2120e13374d43d8a1eba7c46a78da1371234a2f3..55fc32038ba8d51ed7c9c598b051fb16ca22e5e8 100644
--- a/MaterialLib/CMakeLists.txt
+++ b/MaterialLib/CMakeLists.txt
@@ -2,6 +2,7 @@
get_source_files(SOURCES)
append_source_files(SOURCES Adsorption)
append_source_files(SOURCES SolidModels)
+append_source_files(SOURCES FractureModels)
append_source_files(SOURCES Fluid)
append_source_files(SOURCES Fluid/Density)
diff --git a/MaterialLib/FractureModels/CreateLinearElasticIsotropic.cpp b/MaterialLib/FractureModels/CreateLinearElasticIsotropic.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7ae2292d0a88cfb0452eed0cf75014431d811ddc
--- /dev/null
+++ b/MaterialLib/FractureModels/CreateLinearElasticIsotropic.cpp
@@ -0,0 +1,56 @@
+/**
+ * \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 "CreateLinearElasticIsotropic.h"
+
+#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
+#include "LinearElasticIsotropic.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+std::unique_ptr<FractureModelBase<DisplacementDim>>
+createLinearElasticIsotropic(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config)
+{
+ config.checkConfigParameter("type", "LinearElasticIsotropic");
+ DBUG("Create LinearElasticIsotropic material");
+
+ auto& Kn = ProcessLib::findParameter<double>(
+ config, "normal_stiffness", parameters, 1);
+
+ auto& Ks = ProcessLib::findParameter<double>(
+ config, "shear_stiffness", parameters, 1);
+
+ typename LinearElasticIsotropic<DisplacementDim>::MaterialProperties mp{
+ Kn, Ks};
+
+ return std::unique_ptr<LinearElasticIsotropic<DisplacementDim>>{
+ new LinearElasticIsotropic<DisplacementDim>{mp}};
+}
+
+
+template
+std::unique_ptr<FractureModelBase<2>>
+createLinearElasticIsotropic(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+template
+std::unique_ptr<FractureModelBase<3>>
+createLinearElasticIsotropic(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+} // namespace Fracture
+} // namespace MaterialLib
diff --git a/MaterialLib/FractureModels/CreateLinearElasticIsotropic.h b/MaterialLib/FractureModels/CreateLinearElasticIsotropic.h
new file mode 100644
index 0000000000000000000000000000000000000000..1b2e0c4d7ec5e2473fbd33d6e27cd568b7571c9f
--- /dev/null
+++ b/MaterialLib/FractureModels/CreateLinearElasticIsotropic.h
@@ -0,0 +1,29 @@
+/**
+ * \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_FRACTURE_CREATELINEARELASTICISOTROPIC_H_
+#define MATERIALLIB_FRACTURE_CREATELINEARELASTICISOTROPIC_H_
+
+#include "FractureModelBase.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+std::unique_ptr<FractureModelBase<DisplacementDim>>
+createLinearElasticIsotropic(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+} // namespace Fracture
+} // namespace MaterialLib
+
+#endif
diff --git a/MaterialLib/FractureModels/CreateMohrCoulomb.cpp b/MaterialLib/FractureModels/CreateMohrCoulomb.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8e9e546389e80dc410176839d9e2536a3fbece33
--- /dev/null
+++ b/MaterialLib/FractureModels/CreateMohrCoulomb.cpp
@@ -0,0 +1,66 @@
+/**
+ * \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 "CreateMohrCoulomb.h"
+
+#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
+#include "MohrCoulomb.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+std::unique_ptr<FractureModelBase<DisplacementDim>>
+createMohrCoulomb(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config)
+{
+ config.checkConfigParameter("type", "MohrCoulomb");
+ DBUG("Create MohrCoulomb material");
+
+ auto& Kn = ProcessLib::findParameter<double>(
+ config, "normal_stiffness", parameters, 1);
+
+ auto& Ks = ProcessLib::findParameter<double>(
+ config, "shear_stiffness", parameters, 1);
+
+ auto& friction_angle = ProcessLib::findParameter<double>(
+ config, "friction_angle", parameters, 1);
+
+ auto& dilatancy_angle = ProcessLib::findParameter<double>(
+ config, "dilatancy_angle", parameters, 1);
+
+ auto& cohesion = ProcessLib::findParameter<double>(
+ config, "cohesion", parameters, 1);
+
+
+ typename MohrCoulomb<DisplacementDim>::MaterialProperties mp{
+ Kn, Ks, friction_angle, dilatancy_angle, cohesion};
+
+ return std::unique_ptr<MohrCoulomb<DisplacementDim>>{
+ new MohrCoulomb<DisplacementDim>{mp}};
+}
+
+
+template
+std::unique_ptr<FractureModelBase<2>>
+createMohrCoulomb(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+template
+std::unique_ptr<FractureModelBase<3>>
+createMohrCoulomb(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+} // namespace Fracture
+} // namespace MaterialLib
diff --git a/MaterialLib/FractureModels/CreateMohrCoulomb.h b/MaterialLib/FractureModels/CreateMohrCoulomb.h
new file mode 100644
index 0000000000000000000000000000000000000000..1d05598760c681bb717c6d83975b04cf847ca867
--- /dev/null
+++ b/MaterialLib/FractureModels/CreateMohrCoulomb.h
@@ -0,0 +1,29 @@
+/**
+ * \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_FRACTURE_CREATEMOHRCOULOMB_H_
+#define MATERIALLIB_FRACTURE_CREATEMOHRCOULOMB_H_
+
+#include "FractureModelBase.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+std::unique_ptr<FractureModelBase<DisplacementDim>>
+createMohrCoulomb(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+} // namespace Fracture
+} // namespace MaterialLib
+
+#endif
diff --git a/MaterialLib/FractureModels/FractureModelBase.h b/MaterialLib/FractureModels/FractureModelBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..7bcd02e4a640f077eb14f9c4b0f796cd2a698fb8
--- /dev/null
+++ b/MaterialLib/FractureModels/FractureModelBase.h
@@ -0,0 +1,58 @@
+/**
+ * \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_FRACTURE_FRACTUREMODELBASE_H_
+#define MATERIALLIB_FRACTURE_FRACTUREMODELBASE_H_
+
+#include <Eigen/Eigen>
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+/**
+ * Interface for mechanical fracture material models. Provides updates of the
+ * stress for a given current state and also a tangent at that position.
+ */
+template <int DisplacementDim>
+class FractureModelBase
+{
+public:
+ virtual ~FractureModelBase() {}
+
+ /**
+ * Computation of the constitutive relation for specific material model.
+ * This should be implemented in the derived model.
+ *
+ * @param t current time
+ * @param x current position in space
+ * @param w_prev fracture displacement at previous time step
+ * @param w fracture displacement at current time step
+ * @param sigma_prev stress at previous time step
+ * @param sigma stress at current time step
+ * @param C tangent matrix for stress and fracture displacements
+ */
+ virtual void computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ Eigen::Ref<Eigen::VectorXd const> w_prev,
+ Eigen::Ref<Eigen::VectorXd const> w,
+ Eigen::Ref<Eigen::VectorXd const> sigma_prev,
+ Eigen::Ref<Eigen::VectorXd> sigma,
+ Eigen::Ref<Eigen::MatrixXd> C) = 0;
+
+};
+
+} // namespace Fracture
+} // namespace MaterialLib
+
+#endif
diff --git a/MaterialLib/FractureModels/LinearElasticIsotropic.cpp b/MaterialLib/FractureModels/LinearElasticIsotropic.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b7fbf541e0382c080d7cfae6676bfb71a8032a42
--- /dev/null
+++ b/MaterialLib/FractureModels/LinearElasticIsotropic.cpp
@@ -0,0 +1,48 @@
+/**
+ * \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 "LinearElasticIsotropic.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+void LinearElasticIsotropic<DisplacementDim>::computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ Eigen::Ref<Eigen::VectorXd const> w_prev,
+ Eigen::Ref<Eigen::VectorXd const> w,
+ Eigen::Ref<Eigen::VectorXd const> sigma_prev,
+ Eigen::Ref<Eigen::VectorXd> sigma,
+ Eigen::Ref<Eigen::MatrixXd> C)
+{
+ const int index_ns = DisplacementDim - 1;
+ C.setZero();
+ for (int i=0; i<index_ns; i++)
+ C(i,i) = _mp.shear_stiffness(t, x)[0];
+ C(index_ns, index_ns) = _mp.normal_stiffness(t, x)[0];
+
+ sigma.noalias() = sigma_prev + C * (w - w_prev);
+
+ // correct if fracture is opening
+ if (sigma[index_ns] > 0)
+ {
+ C.setZero();
+ sigma.setZero();
+ }
+}
+
+template class LinearElasticIsotropic<2>;
+template class LinearElasticIsotropic<3>;
+
+
+} // namespace Fracture
+} // namespace MaterialLib
diff --git a/MaterialLib/FractureModels/LinearElasticIsotropic.h b/MaterialLib/FractureModels/LinearElasticIsotropic.h
new file mode 100644
index 0000000000000000000000000000000000000000..c236be4e10b34a45f10d0b99867376d7d31a19fc
--- /dev/null
+++ b/MaterialLib/FractureModels/LinearElasticIsotropic.h
@@ -0,0 +1,75 @@
+/**
+ * \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_FRACTURE_LINEARELASTICISOTROPIC_H_
+#define MATERIALLIB_FRACTURE_LINEARELASTICISOTROPIC_H_
+
+#include <Eigen/Eigen>
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+#include "FractureModelBase.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+class LinearElasticIsotropic final : public FractureModelBase<DisplacementDim>
+{
+public:
+ /// Variables specific to the material model
+ struct MaterialProperties
+ {
+ using P = ProcessLib::Parameter<double>;
+ using X = ProcessLib::SpatialPosition;
+
+ MaterialProperties(P const& normal_stiffness_, P const& shear_stiffness_)
+ : normal_stiffness(normal_stiffness_), shear_stiffness(shear_stiffness_)
+ {
+ }
+
+ P const& normal_stiffness;
+ P const& shear_stiffness;
+ };
+
+public:
+ explicit LinearElasticIsotropic(
+ MaterialProperties const& material_properties)
+ : _mp(material_properties)
+ {
+ }
+
+ void computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ Eigen::Ref<Eigen::VectorXd const> w_prev,
+ Eigen::Ref<Eigen::VectorXd const> w,
+ Eigen::Ref<Eigen::VectorXd const> sigma_prev,
+ Eigen::Ref<Eigen::VectorXd> sigma,
+ Eigen::Ref<Eigen::MatrixXd> C) override;
+
+private:
+ MaterialProperties _mp;
+};
+
+} // namespace Fracture
+} // namespace MaterialLib
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+extern template class LinearElasticIsotropic<2>;
+extern template class LinearElasticIsotropic<3>;
+} // namespace Fractrue
+} // namespace MaterialLib
+
+#endif
diff --git a/MaterialLib/FractureModels/MohrCoulomb.cpp b/MaterialLib/FractureModels/MohrCoulomb.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..175be798f772e0654678c2d573d0c95a2b946f94
--- /dev/null
+++ b/MaterialLib/FractureModels/MohrCoulomb.cpp
@@ -0,0 +1,119 @@
+/**
+ * \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 "MohrCoulomb.h"
+
+#include "BaseLib/Error.h"
+#include "MathLib/MathTools.h"
+
+#ifndef Q_MOC_RUN // to avoid Qt4 bug, https://bugreports.qt.io/browse/QTBUG-22829
+#include <boost/math/special_functions/sign.hpp>
+#endif
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+namespace
+{
+
+struct MaterialPropertyValues
+{
+ double Kn = 0.0;
+ double Ks = 0.0;
+ double phi = 0.0; // friction angle
+ double psi = 0.0; // dilation angle
+ double c = 0.0;
+
+ template <typename MaterialProperties>
+ MaterialPropertyValues(
+ MaterialProperties const& mp,
+ double const t,
+ ProcessLib::SpatialPosition const& x)
+ {
+ Kn = mp.normal_stiffness(t,x)[0];
+ Ks = mp.shear_stiffness(t,x)[0];
+ phi = MathLib::to_radians(mp.friction_angle(t,x)[0]);
+ psi = MathLib::to_radians(mp.dilatancy_angle(t,x)[0]);
+ c = mp.cohesion(t,x)[0];
+ }
+};
+
+} // no namespace
+
+template <int DisplacementDim>
+void MohrCoulomb<DisplacementDim>::computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ Eigen::Ref<Eigen::VectorXd const> w_prev,
+ Eigen::Ref<Eigen::VectorXd const> w,
+ Eigen::Ref<Eigen::VectorXd const> sigma_prev,
+ Eigen::Ref<Eigen::VectorXd> sigma,
+ Eigen::Ref<Eigen::MatrixXd> Kep)
+{
+ if (DisplacementDim == 3)
+ {
+ OGS_FATAL("MohrCoulomb fracture model does not support 3D case.");
+ return;
+ }
+ MaterialPropertyValues const mat(_mp, t, x);
+ Eigen::VectorXd const dw = w - w_prev;
+
+ Eigen::MatrixXd Ke(2,2);
+ Ke.setZero();
+ Ke(0,0) = mat.Ks;
+ Ke(1,1) = mat.Kn;
+
+ sigma.noalias() = sigma_prev + Ke * dw;
+
+ // if opening
+ if (sigma[1] > 0)
+ {
+ Kep.setZero();
+ sigma.setZero();
+ return;
+ }
+
+ // check shear yield function (Fs)
+ double const Fs = std::abs(sigma[0]) + sigma[1] * std::tan(mat.phi) - mat.c;
+ if (Fs < .0)
+ {
+ Kep = Ke;
+ return;
+ }
+
+ Eigen::VectorXd dFs_dS(2);
+ dFs_dS[0] = boost::math::sign(sigma[0]);
+ dFs_dS[1] = std::tan(mat.phi);
+
+ // plastic potential function: Qs = |tau| + Sn * tan da
+ Eigen::VectorXd dQs_dS(2);
+ dQs_dS[0] = boost::math::sign(sigma[0]);
+ dQs_dS[1] = std::tan(mat.psi);
+
+ // plastic multiplier
+ Eigen::RowVectorXd const A = dFs_dS.transpose() * Ke / (dFs_dS.transpose() * Ke * dQs_dS);
+ double const d_eta = A * dw;
+
+ // plastic part of the dispalcement
+ Eigen::VectorXd const dwp = dQs_dS * d_eta;
+
+ // correct stress
+ sigma.noalias() = sigma_prev + Ke * (dw - dwp);
+
+ // Kep
+ Kep = Ke - Ke * dQs_dS * A;
+}
+
+template class MohrCoulomb<2>;
+template class MohrCoulomb<3>;
+
+} // namespace Fracture
+} // namespace MaterialLib
diff --git a/MaterialLib/FractureModels/MohrCoulomb.h b/MaterialLib/FractureModels/MohrCoulomb.h
new file mode 100644
index 0000000000000000000000000000000000000000..cd7d974ace555ba92bffcbbf40060ab4c6b9caed
--- /dev/null
+++ b/MaterialLib/FractureModels/MohrCoulomb.h
@@ -0,0 +1,85 @@
+/**
+ * \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_FRACTURE_MOHRCOULOMB_H_
+#define MATERIALLIB_FRACTURE_MOHRCOULOMB_H_
+
+#include <Eigen/Eigen>
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+#include "FractureModelBase.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+class MohrCoulomb final : public FractureModelBase<DisplacementDim>
+{
+public:
+ /// Variables specific to the material model
+ struct MaterialProperties
+ {
+ using P = ProcessLib::Parameter<double>;
+ using X = ProcessLib::SpatialPosition;
+
+ MaterialProperties(
+ P const& normal_stiffness_, P const& shear_stiffness_,
+ P const& friction_angle_, P const& dilatancy_angle_,
+ P const& cohesion_)
+ : normal_stiffness(normal_stiffness_), shear_stiffness(shear_stiffness_),
+ friction_angle(friction_angle_), dilatancy_angle(dilatancy_angle_),
+ cohesion(cohesion_)
+ {
+ }
+
+ P const& normal_stiffness;
+ P const& shear_stiffness;
+ P const& friction_angle;
+ P const& dilatancy_angle;
+ P const& cohesion;
+ };
+
+public:
+
+ explicit MohrCoulomb(
+ MaterialProperties const& material_properties)
+ : _mp(material_properties)
+ {
+ }
+
+ void computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ Eigen::Ref<Eigen::VectorXd const> w_prev,
+ Eigen::Ref<Eigen::VectorXd const> w,
+ Eigen::Ref<Eigen::VectorXd const> sigma_prev,
+ Eigen::Ref<Eigen::VectorXd> sigma,
+ Eigen::Ref<Eigen::MatrixXd> Kep) override;
+
+private:
+
+ MaterialProperties _mp;
+};
+
+} // namespace Fracture
+} // namespace MaterialLib
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+extern template class MohrCoulomb<2>;
+extern template class MohrCoulomb<3>;
+} // namespace Fractrue
+} // namespace MaterialLib
+
+#endif
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index 81f1fad011656ef37d106cf515fe5b4c6960c99c..e578be54451e0447bbd9aa198bf93c8c2c21140e 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -17,6 +17,10 @@
#include <omp.h>
#endif
+#ifndef Q_MOC_RUN // to avoid Qt4 bug, https://bugreports.qt.io/browse/QTBUG-22829
+#include <boost/math/constants/constants.hpp>
+#endif
+
namespace MathLib
{
/**
@@ -111,6 +115,11 @@ double sqrDist(const double* p0, const double* p1)
*/
double getAngle (const double p0[3], const double p1[3], const double p2[3]);
+/// converts the given degrees to radians
+inline double to_radians(double degrees) {
+ return degrees*boost::math::constants::pi<double>()/180.;
+}
+
} // namespace
#endif /* MATHTOOLS_H_ */
diff --git a/Tests/MaterialLib/TestFractureModels.cpp b/Tests/MaterialLib/TestFractureModels.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..18857a286c060bddcc99049e33e99e2d63044c38
--- /dev/null
+++ b/Tests/MaterialLib/TestFractureModels.cpp
@@ -0,0 +1,88 @@
+/**
+ * \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/LICENSE.txt
+ */
+
+
+#include <limits>
+
+#include <gtest/gtest.h>
+
+#include <Eigen/Eigen>
+
+#include "BaseLib/ConfigTree.h"
+
+#include "MaterialLib/FractureModels/LinearElasticIsotropic.h"
+#include "MaterialLib/FractureModels/MohrCoulomb.h"
+
+#include "ProcessLib/Parameter/ConstantParameter.h"
+
+using namespace MaterialLib::Fracture;
+
+static const double eps_sigma = 1e6*1e-5;
+static const double eps_C = 1e10*1e-5;
+
+TEST(MaterialLib_Fracture, LinearElasticIsotropic)
+{
+ ProcessLib::ConstantParameter<double> const kn(1e11);
+ ProcessLib::ConstantParameter<double> const ks(1e9);
+ LinearElasticIsotropic<2>::MaterialProperties const mp{kn, ks};
+
+ LinearElasticIsotropic<2> fractureModel{mp};
+
+ Eigen::Vector2d w_prev, w, sigma_prev, sigma;
+ Eigen::Matrix2d C;
+
+ ProcessLib::SpatialPosition x;
+ w << -1e-5, -1e-5;
+ fractureModel.computeConstitutiveRelation(0, x, w_prev, w, sigma_prev, sigma, C);
+
+ ASSERT_NEAR(-1e4, sigma[0], eps_sigma);
+ ASSERT_NEAR(-1e6, sigma[1], eps_sigma);
+ ASSERT_NEAR(1e9, C(0,0), eps_C);
+ ASSERT_NEAR(0, C(0,1), eps_C);
+ ASSERT_NEAR(0, C(1,0), eps_C);
+ ASSERT_NEAR(1e11, C(1,1), eps_C);
+
+
+ w << -1e-5, 1e-5;
+ fractureModel.computeConstitutiveRelation(0, x, w_prev, w, sigma_prev, sigma, C);
+
+ ASSERT_NEAR(0, sigma[0], eps_sigma);
+ ASSERT_NEAR(0, sigma[1], eps_sigma);
+ ASSERT_NEAR(0, C(0,0), eps_C);
+ ASSERT_NEAR(0, C(0,1), eps_C);
+ ASSERT_NEAR(0, C(1,0), eps_C);
+ ASSERT_NEAR(0, C(1,1), eps_C);
+}
+
+TEST(MaterialLib_Fracture, MohrCoulomb)
+{
+ ProcessLib::ConstantParameter<double> const kn(50e9);
+ ProcessLib::ConstantParameter<double> const ks(20e9);
+ ProcessLib::ConstantParameter<double> const phi(15);
+ ProcessLib::ConstantParameter<double> const psi(5);
+ ProcessLib::ConstantParameter<double> const c(3e6);
+ MohrCoulomb<2>::MaterialProperties const mp{kn, ks, phi, psi, c};
+
+ MohrCoulomb<2> fractureModel{mp};
+
+ Eigen::Vector2d w_prev, w, sigma_prev, sigma;
+ Eigen::Matrix2d C;
+
+ ProcessLib::SpatialPosition x;
+ sigma_prev << -3.46e6, -2e6;
+ w << -1.08e-5, -0.25e-5;
+ fractureModel.computeConstitutiveRelation(0, x, w_prev, w, sigma_prev, sigma, C);
+
+ ASSERT_NEAR(-3.50360e6, sigma[0], eps_sigma);
+ ASSERT_NEAR(-2.16271e6, sigma[1], eps_sigma);
+ ASSERT_NEAR(1.10723e+09, C(0,0), eps_C);
+ ASSERT_NEAR(1.26558e+10, C(0,1), eps_C);
+ ASSERT_NEAR(4.13226e+09, C(1,0), eps_C);
+ ASSERT_NEAR(4.72319e+10, C(1,1), eps_C);
+}
+