diff --git a/Documentation/ProjectFile/material/fracture_model/MohrCoulomb/nonlinear_solver b/Documentation/ProjectFile/material/fracture_model/MohrCoulomb/nonlinear_solver
new file mode 120000
index 0000000000000000000000000000000000000000..d410878dc905bfe81ce70c44a5e4e19101394641
--- /dev/null
+++ b/Documentation/ProjectFile/material/fracture_model/MohrCoulomb/nonlinear_solver
@@ -0,0 +1 @@
+../../../nonlinear_solver
\ No newline at end of file
diff --git a/MaterialLib/FractureModels/CreateMohrCoulomb.cpp b/MaterialLib/FractureModels/CreateMohrCoulomb.cpp
index 0f27638fd98e3847d8c075fc392976426cca415b..874695927ed057df2f6c9f04eda318fa8c667e3d 100644
--- a/MaterialLib/FractureModels/CreateMohrCoulomb.cpp
+++ b/MaterialLib/FractureModels/CreateMohrCoulomb.cpp
@@ -9,6 +9,7 @@
#include "CreateMohrCoulomb.h"
+#include "MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.h"
#include "ParameterLib/Utils.h"
#include "MohrCoulomb.h"
@@ -57,8 +58,16 @@ std::unique_ptr<FractureModelBase<DisplacementDim>> createMohrCoulomb(
typename MohrCoulomb::MohrCoulomb<DisplacementDim>::MaterialProperties mp{
Kn, Ks, friction_angle, dilatancy_angle, cohesion};
+ auto const& nonlinear_solver_config =
+ //! \ogs_file_param{material__fracture_model__MohrCoulomb__nonlinear_solver}
+ config.getConfigSubtree("nonlinear_solver");
+ auto const nonlinear_solver_parameters =
+ MaterialLib::createNewtonRaphsonSolverParameters(
+ nonlinear_solver_config);
+
return std::make_unique<MohrCoulomb::MohrCoulomb<DisplacementDim>>(
- penalty_aperture_cutoff, tension_cutoff, mp);
+ nonlinear_solver_parameters, penalty_aperture_cutoff, tension_cutoff,
+ mp);
}
template std::unique_ptr<FractureModelBase<2>> createMohrCoulomb(
diff --git a/MaterialLib/FractureModels/MohrCoulomb.cpp b/MaterialLib/FractureModels/MohrCoulomb.cpp
index cd5ca1abb9eb965a271241103933138eac13faf9..2e7511ff40665978f7774919e6f3ca4b9496fe5a 100644
--- a/MaterialLib/FractureModels/MohrCoulomb.cpp
+++ b/MaterialLib/FractureModels/MohrCoulomb.cpp
@@ -53,11 +53,11 @@ void MohrCoulomb<DisplacementDim>::computeConstitutiveRelation(
Eigen::Ref<Eigen::VectorXd const>
sigma0,
Eigen::Ref<Eigen::VectorXd const>
- w_prev,
+ /*w_prev*/,
Eigen::Ref<Eigen::VectorXd const>
w,
Eigen::Ref<Eigen::VectorXd const>
- sigma_prev,
+ /*sigma_prev*/,
Eigen::Ref<Eigen::VectorXd>
sigma,
Eigen::Ref<Eigen::MatrixXd>
@@ -65,14 +65,16 @@ void MohrCoulomb<DisplacementDim>::computeConstitutiveRelation(
typename FractureModelBase<DisplacementDim>::MaterialStateVariables&
material_state_variables)
{
- material_state_variables.reset();
+ assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
+ &material_state_variables) != nullptr);
+
+ StateVariables<DisplacementDim>& state =
+ static_cast<StateVariables<DisplacementDim>&>(material_state_variables);
MaterialPropertyValues const mat(_mp, t, x);
- Eigen::VectorXd const dw = w - w_prev;
const int index_ns = DisplacementDim - 1;
double const aperture = w[index_ns] + aperture0;
- double const aperture_prev = w_prev[index_ns] + aperture0;
Eigen::MatrixXd Ke;
{ // Elastic tangent stiffness
@@ -82,84 +84,132 @@ void MohrCoulomb<DisplacementDim>::computeConstitutiveRelation(
Ke(i, i) = mat.Ks;
}
- Ke(index_ns, index_ns) =
- mat.Kn *
- logPenaltyDerivative(aperture0, aperture, _penalty_aperture_cutoff);
- }
-
- Eigen::MatrixXd Ke_prev;
- { // Elastic tangent stiffness at w_prev
- Ke_prev = Eigen::MatrixXd::Zero(DisplacementDim, DisplacementDim);
- for (int i = 0; i < index_ns; i++)
- {
- Ke_prev(i, i) = mat.Ks;
- }
-
- Ke_prev(index_ns, index_ns) =
- mat.Kn * logPenaltyDerivative(
- aperture0, aperture_prev, _penalty_aperture_cutoff);
+ Ke(index_ns, index_ns) = mat.Kn;
}
// Total plastic aperture compression
// NOTE: Initial condition sigma0 seems to be associated with an initial
// condition of the w0 = 0. Therefore the initial state is not associated
// with a plastic aperture change.
- Eigen::VectorXd const w_p_prev =
- w_prev - Ke_prev.fullPivLu().solve(sigma_prev - sigma0);
-
{ // Exact elastic predictor
- sigma.noalias() = Ke * (w - w_p_prev);
+ sigma.noalias() = Ke * (w - state.w_p_prev);
- sigma.coeffRef(index_ns) =
- mat.Kn * w[index_ns] *
+ sigma.coeffRef(index_ns) *=
logPenalty(aperture0, aperture, _penalty_aperture_cutoff);
}
sigma.noalias() += sigma0;
- double const sigma_n = sigma[index_ns];
-
// correction for an opening fracture
- if (_tension_cutoff && sigma_n > 0)
+ if (_tension_cutoff && sigma[DisplacementDim - 1] >= 0)
{
Kep.setZero();
sigma.setZero();
material_state_variables.setTensileStress(true);
return;
- }
- // check shear yield function (Fs)
- Eigen::VectorXd const sigma_s = sigma.head(DisplacementDim-1);
- double const mag_tau = sigma_s.norm(); // magnitude
- double const Fs = mag_tau + sigma_n * std::tan(mat.phi) - mat.c;
+ // TODO (nagel); Update w_p for fracture opening and closing.
+ }
- material_state_variables.setShearYieldFunctionValue(Fs);
- if (Fs < .0)
- {
- Kep = Ke;
- return;
+ auto yield_function = [&mat](Eigen::VectorXd const& s) {
+ double const sigma_n = s[DisplacementDim - 1];
+ Eigen::VectorXd const sigma_s = s.head(DisplacementDim - 1);
+ double const mag_tau = sigma_s.norm(); // magnitude
+ return mag_tau + sigma_n * std::tan(mat.phi) - mat.c;
+ };
+
+ { // Exit if still in elastic range by checking the shear yield function.
+ double const Fs = yield_function(sigma);
+ material_state_variables.setShearYieldFunctionValue(Fs);
+ if (Fs < .0)
+ {
+ Kep = Ke;
+ Kep(index_ns, index_ns) *= logPenaltyDerivative(
+ aperture0, aperture, _penalty_aperture_cutoff);
+ return;
+ }
}
- Eigen::VectorXd dFs_dS(DisplacementDim);
- dFs_dS.head(DisplacementDim-1).noalias() = sigma_s.normalized();
- dFs_dS[index_ns] = std::tan(mat.phi);
+ auto yield_function_derivative = [&mat](Eigen::VectorXd const& s) {
+ Eigen::Matrix<double, DisplacementDim, 1> dFs_dS;
+ dFs_dS.template head<DisplacementDim - 1>().noalias() =
+ s.template head<DisplacementDim - 1>().normalized();
+ dFs_dS.coeffRef(DisplacementDim - 1) = std::tan(mat.phi);
+ return dFs_dS;
+ };
// plastic potential function: Qs = |tau| + Sn * tan da
- Eigen::VectorXd dQs_dS = dFs_dS;
- dQs_dS[index_ns] = 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;
+ auto plastic_potential_derivative = [&mat](Eigen::VectorXd const& s) {
+ Eigen::Matrix<double, DisplacementDim, 1> dQs_dS;
+ dQs_dS.template head<DisplacementDim - 1>().noalias() =
+ s.template head<DisplacementDim - 1>().normalized();
+ dQs_dS.coeffRef(DisplacementDim - 1) = std::tan(mat.psi);
+ return dQs_dS;
+ };
+
+ { // Newton
+
+ Eigen::FullPivLU<Eigen::Matrix<double, 1, 1, Eigen::RowMajor>>
+ linear_solver;
+ using ResidualVectorType = Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
+ using JacobianMatrix = Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
+
+ JacobianMatrix jacobian;
+ ResidualVectorType solution;
+ solution << 0;
+
+ auto const update_residual = [&](ResidualVectorType& residual) {
+ residual[0] = yield_function(sigma);
+ };
+
+ auto const update_jacobian = [&](JacobianMatrix& jacobian) {
+ jacobian(0, 0) = -yield_function_derivative(sigma).transpose() *
+ Ke * plastic_potential_derivative(sigma);
+ };
+
+ auto const update_solution = [&](ResidualVectorType const& increment) {
+ solution += increment;
+ /*DBUG("analytical = %g",
+ Fs / (mat.Ks + mat.Kn * std::tan(mat.psi) * std::tan(mat.phi)))
+ */
+ state.w_p = state.w_p_prev +
+ solution[0] * plastic_potential_derivative(sigma);
+ sigma.noalias() = sigma0 + (Ke * (w - state.w_p)) *
+ logPenalty(aperture0, aperture,
+ _penalty_aperture_cutoff);
+ };
+
+ auto newton_solver =
+ NumLib::NewtonRaphson<decltype(linear_solver), JacobianMatrix,
+ decltype(update_jacobian), ResidualVectorType,
+ decltype(update_residual),
+ decltype(update_solution)>(
+ linear_solver, update_jacobian, update_residual,
+ update_solution, _nonlinear_solver_parameters);
+
+ auto const success_iterations = newton_solver.solve(jacobian);
+
+ if (!success_iterations)
+ {
+ OGS_FATAL("MohrCoulomb nonlinear solver didn't converge.");
+ }
- // plastic part of the dispalcement
- Eigen::VectorXd const dwp = dQs_dS * d_eta;
+ // Solution containing lambda is not needed; w_p and sigma already
+ // up to date.
+ }
- // correct stress
- sigma.noalias() = sigma_prev + Ke * (dw - dwp);
+ { // Update material state shear yield function value.
+ double const Fs = yield_function(sigma);
+ material_state_variables.setShearYieldFunctionValue(Fs);
+ }
- // Kep
- Kep = Ke - Ke * dQs_dS * A;
+ Ke(index_ns, index_ns) *=
+ logPenaltyDerivative(aperture0, aperture, _penalty_aperture_cutoff);
+ Eigen::RowVectorXd const A = yield_function_derivative(sigma).transpose() *
+ Ke /
+ (yield_function_derivative(sigma).transpose() *
+ Ke * plastic_potential_derivative(sigma));
+ Kep = Ke - Ke * plastic_potential_derivative(sigma) * A;
}
template class MohrCoulomb<2>;
diff --git a/MaterialLib/FractureModels/MohrCoulomb.h b/MaterialLib/FractureModels/MohrCoulomb.h
index 3b9e522e285d6eff975277b67e533e51f1f69689..c0e2ed2234a9a55215808fd7f47257d2834e8f5c 100644
--- a/MaterialLib/FractureModels/MohrCoulomb.h
+++ b/MaterialLib/FractureModels/MohrCoulomb.h
@@ -12,6 +12,7 @@
#include <Eigen/Eigen>
#include <utility>
+#include "NumLib/NewtonRaphson.h"
#include "ParameterLib/Parameter.h"
#include "FractureModelBase.h"
@@ -22,9 +23,36 @@ namespace Fracture
{
namespace MohrCoulomb
{
+template <int DisplacementDim>
+struct StateVariables
+ : public FractureModelBase<DisplacementDim>::MaterialStateVariables
+{
+ void setInitialConditions() { w_p = w_p_prev; }
+
+ void pushBackState() override { w_p_prev = w_p; }
+
+ /// Plastic component of the displacement jump in fracture's local
+ /// coordinates.
+ Eigen::Matrix<double, DisplacementDim, 1> w_p =
+ Eigen::Matrix<double, DisplacementDim, 1>::Zero();
+
+ // Initial values from previous timestep
+ Eigen::Matrix<double, DisplacementDim, 1> w_p_prev; ///< \copydoc w_p
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
template <int DisplacementDim>
class MohrCoulomb final : public FractureModelBase<DisplacementDim>
{
+public:
+ std::unique_ptr<
+ typename FractureModelBase<DisplacementDim>::MaterialStateVariables>
+ createMaterialStateVariables() override
+ {
+ return std::make_unique<StateVariables<DisplacementDim>>();
+ }
+
public:
/// Variables specific to the material model
struct MaterialProperties
@@ -58,26 +86,14 @@ public:
};
- struct MaterialStateVariables
- : public FractureModelBase<DisplacementDim>::MaterialStateVariables
- {
- void pushBackState() override {}
- };
-
- std::unique_ptr<
- typename FractureModelBase<DisplacementDim>::MaterialStateVariables>
- createMaterialStateVariables() override
- {
- return std::unique_ptr<typename FractureModelBase<
- DisplacementDim>::MaterialStateVariables>{
- new MaterialStateVariables};
- }
-
public:
- explicit MohrCoulomb(double const penalty_aperture_cutoff,
- bool const tension_cutoff,
- MaterialProperties material_properties)
- : _penalty_aperture_cutoff(penalty_aperture_cutoff),
+ explicit MohrCoulomb(
+ NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters,
+ double const penalty_aperture_cutoff,
+ bool const tension_cutoff,
+ MaterialProperties material_properties)
+ : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
+ _penalty_aperture_cutoff(penalty_aperture_cutoff),
_tension_cutoff(tension_cutoff),
_mp(std::move(material_properties))
{
@@ -117,6 +133,8 @@ public:
material_state_variables) override;
private:
+ NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters;
+
/// Compressive normal displacements above this value will not enter the
/// computation of the normal stiffness modulus of the fracture.
/// \note Setting this to the initial aperture value allows negative
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture.h
index f473ad8e3a33e2b37b114ea9616e4a12919d3e32..be07a6d3fbd7eb0179d5fcda6fcc9c08b4aed364 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrixNearFracture.h
@@ -60,6 +60,16 @@ private:
Eigen::VectorXd& local_b,
Eigen::MatrixXd& local_J) override;
+ void preTimestepConcrete(std::vector<double> const& /*local_x*/,
+ double const /*t*/,
+ double const /*delta_t*/) override
+ {
+ for (auto& ip_data : _ip_data)
+ {
+ ip_data.pushBackState();
+ }
+ }
+
void computeSecondaryVariableConcreteWithVector(
double const t, Eigen::VectorXd const& local_x) override;
diff --git a/ProcessLib/LIE/SmallDeformation/Tests.cmake b/ProcessLib/LIE/SmallDeformation/Tests.cmake
index 39b380a42ee9102688c14a16dc16b91f47f98e0f..0dffe32ccd639db97c669f88238eeb3d1440f0f0 100644
--- a/ProcessLib/LIE/SmallDeformation/Tests.cmake
+++ b/ProcessLib/LIE/SmallDeformation/Tests.cmake
@@ -75,8 +75,8 @@ AddTest(
TESTER vtkdiff
REQUIREMENTS NOT OGS_USE_MPI
DIFF_DATA
- expected_single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu displacement displacement 1e-16 1e-16
- expected_single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu displacement_jump1 displacement_jump1 1e-16 1e-16
+ expected_single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu displacement displacement 1e-16 0
+ expected_single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu displacement_jump1 displacement_jump1 1e-16 0
expected_single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu aperture aperture 1e-16 1e-16
)
diff --git a/Tests/Data/LIE/Mechanics/expected_mohr_coulomb_load_path_pcs_0_ts_100_t_1.000000.vtu b/Tests/Data/LIE/Mechanics/expected_mohr_coulomb_load_path_pcs_0_ts_100_t_1.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..9bf8c80855eb263e9dff9d4be37e4b00403d7d47
--- /dev/null
+++ b/Tests/Data/LIE/Mechanics/expected_mohr_coulomb_load_path_pcs_0_ts_100_t_1.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:71a3b44705cc0e527e05526b352edd1c6d4918e2b4ae6a3a505ce40655264287
+size 4433
diff --git a/Tests/Data/LIE/Mechanics/expected_mohr_coulomb_load_path_pcs_0_ts_200_t_2.000000.vtu b/Tests/Data/LIE/Mechanics/expected_mohr_coulomb_load_path_pcs_0_ts_200_t_2.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..90c878be402226080855086778fdef7156317daa
--- /dev/null
+++ b/Tests/Data/LIE/Mechanics/expected_mohr_coulomb_load_path_pcs_0_ts_200_t_2.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6d99bd3d2b430b776f420c3700cf4a2080070ba6318f86d4b4f3071716d7bb48
+size 4460
diff --git a/Tests/Data/LIE/Mechanics/expected_mohr_coulomb_load_path_pcs_0_ts_401_t_4.000000.vtu b/Tests/Data/LIE/Mechanics/expected_mohr_coulomb_load_path_pcs_0_ts_401_t_4.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..48d5eac48dbb24fd64b38064048cb691cd2b2587
--- /dev/null
+++ b/Tests/Data/LIE/Mechanics/expected_mohr_coulomb_load_path_pcs_0_ts_401_t_4.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1e265fc5987b57e3236f6a17cd95e273b563bb7e5d499c95ce562eaec769cdee
+size 4478
diff --git a/Tests/Data/LIE/Mechanics/expected_single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu b/Tests/Data/LIE/Mechanics/expected_single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu
index 2582e1f3c71ba4b03dc689213b10963cd2cdba9e..f56fb3a2e9b8ca611e1a6186c43e7bc6f5e1c424 100644
--- a/Tests/Data/LIE/Mechanics/expected_single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu
+++ b/Tests/Data/LIE/Mechanics/expected_single_joint_displacement_controlled_pcs_0_ts_10_t_1.000000.vtu
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:fd3522e936d823a0044c7ac34bab82f31976a7639c5db5cadaf3516871517863
-size 22475
+oid sha256:227d8b102806a342bbf649c94a74475e1305dce8ead099a52a7e911877dfd36b
+size 25478
diff --git a/Tests/Data/LIE/Mechanics/mohr_coulomb_load_path.gml b/Tests/Data/LIE/Mechanics/mohr_coulomb_load_path.gml
new file mode 100644
index 0000000000000000000000000000000000000000..b304ced148e8ab619ffcd159cc710b056eed19b1
--- /dev/null
+++ b/Tests/Data/LIE/Mechanics/mohr_coulomb_load_path.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:05f98a3c2c0fca6c41e900fda62bba5d69562e94e0c0a84ae51639b90a04e2ce
+size 881
diff --git a/Tests/Data/LIE/Mechanics/mohr_coulomb_load_path.prj b/Tests/Data/LIE/Mechanics/mohr_coulomb_load_path.prj
new file mode 100644
index 0000000000000000000000000000000000000000..bc9d676c85df48cd6731ddfaa2290875a6a88270
--- /dev/null
+++ b/Tests/Data/LIE/Mechanics/mohr_coulomb_load_path.prj
@@ -0,0 +1,236 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>mohr_coulomb_load_path.vtu</mesh>
+ <geometry>mohr_coulomb_load_path.gml</geometry>
+ <processes>
+ <process>
+ <name>SD</name>
+ <type>SMALL_DEFORMATION_WITH_LIE</type>
+ <integration_order>2</integration_order>
+ <dimension>2</dimension>
+ <process_variables>
+ <process_variable>displacement</process_variable>
+ <process_variable>displacement_jump1</process_variable>
+ </process_variables>
+ <constitutive_relation>
+ <type>LinearElasticIsotropic</type>
+ <youngs_modulus>E</youngs_modulus>
+ <poissons_ratio>nu</poissons_ratio>
+ </constitutive_relation>
+ <fracture_model>
+ <type>MohrCoulomb</type>
+ <normal_stiffness>Kn</normal_stiffness>
+ <shear_stiffness>Ks</shear_stiffness>
+ <friction_angle>friction_angle</friction_angle>
+ <dilatancy_angle>dilatancy_angle</dilatancy_angle>
+ <cohesion>cohesion</cohesion>
+ <penalty_aperture_cutoff>10e-6</penalty_aperture_cutoff>
+ <tension_cutoff>1</tension_cutoff>
+ <nonlinear_solver>
+ <maximum_iterations>100</maximum_iterations>
+ <error_tolerance>1e-10</error_tolerance>
+ </nonlinear_solver>
+ </fracture_model>
+ <fracture_properties>
+ <material_id>0</material_id>
+ <initial_aperture>aperture0</initial_aperture>
+ </fracture_properties>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="sigma_xx" output_name="sigma_xx"/>
+ <secondary_variable type="static" internal_name="sigma_yy" output_name="sigma_yy"/>
+ <secondary_variable type="static" internal_name="sigma_xy" output_name="sigma_xy"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="SD">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>PerComponentDeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstols>1e-16 1e-16 1e-16 1e-16 1e-16 1e-16</abstols>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>4</t_end>
+ <timesteps>
+ <pair>
+ <repeat>10000</repeat>
+ <delta_t>0.01</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>mohr_coulomb_load_path</prefix>
+ <timesteps>
+ <pair>
+ <repeat>5</repeat>
+ <each_steps>1</each_steps>
+ </pair>
+ </timesteps>
+ <variables>
+ <variable>displacement</variable>
+ <variable>displacement_jump1</variable>
+ <variable>sigma_xx</variable>
+ <variable>sigma_yy</variable>
+ <variable>sigma_xy</variable>
+ </variables>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>E</name>
+ <type>Constant</type>
+ <value>200e+16</value>
+ </parameter>
+ <parameter>
+ <name>nu</name>
+ <type>Constant</type>
+ <value>0.3</value>
+ </parameter>
+ <!-- Fracture properties -->
+ <parameter>
+ <name>Kn</name>
+ <type>Constant</type>
+ <value>20e9</value>
+ </parameter>
+ <parameter>
+ <name>Ks</name>
+ <type>Constant</type>
+ <value>20e9</value>
+ </parameter>
+ <parameter>
+ <name>friction_angle</name>
+ <type>Constant</type>
+ <value>22</value>
+ </parameter>
+ <parameter>
+ <name>dilatancy_angle</name>
+ <type>Constant</type>
+ <value>10</value>
+ </parameter>
+ <parameter>
+ <name>cohesion</name>
+ <type>Constant</type>
+ <value>1e3</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0</values>
+ </parameter>
+ <parameter>
+ <name>aperture0</name>
+ <type>Constant</type>
+ <value>1e-5</value>
+ </parameter>
+ <parameter>
+ <name>displacement_unit</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>displacement_x</name>
+ <type>CurveScaled</type>
+ <curve>timeRamp_displacement_x</curve>
+ <parameter>displacement_unit</parameter>
+ </parameter>
+ <parameter>
+ <name>displacement_y</name>
+ <type>CurveScaled</type>
+ <curve>timeRamp_displacement_y</curve>
+ <parameter>displacement_unit</parameter>
+ </parameter>
+ <parameter>
+ <name>zero_u</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>timeRamp_displacement_x</name>
+ <coords>0 1 2 4 10 </coords>
+ <values>0 0 1e-5 -1e-5 -1e-5</values>
+ </curve>
+ <curve>
+ <name>timeRamp_displacement_y</name>
+ <coords>0 1 10 </coords>
+ <values>0 -1e-6 -1e-6</values>
+ </curve>
+ </curves>
+ <process_variables>
+ <process_variable>
+ <name>displacement</name>
+ <components>2</components>
+ <order>1</order>
+ <initial_condition>displacement0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>DoubleEdgeCrack</geometrical_set>
+ <geometry>PLY_TOP</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>displacement_x</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>DoubleEdgeCrack</geometrical_set>
+ <geometry>PLY_TOP</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>displacement_y</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>DoubleEdgeCrack</geometrical_set>
+ <geometry>PLY_BOTTOM</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>zero_u</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>DoubleEdgeCrack</geometrical_set>
+ <geometry>PLY_BOTTOM</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>zero_u</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>displacement_jump1</name>
+ <components>2</components>
+ <order>1</order>
+ <initial_condition>displacement0</initial_condition>
+ <boundary_conditions>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>25</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <eigen>
+ <solver_type>SparseLU</solver_type>
+ <scaling>1</scaling>
+ <precon_type>DIAGONAL</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/LIE/Mechanics/mohr_coulomb_load_path.vtu b/Tests/Data/LIE/Mechanics/mohr_coulomb_load_path.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..0d37bfd14e7bb42237eca6bdf690ae78c662222d
--- /dev/null
+++ b/Tests/Data/LIE/Mechanics/mohr_coulomb_load_path.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:935a92068d70978a79faaffb8c578e817c55f82bad36df3384ad83c849ca913d
+size 1436
diff --git a/Tests/Data/LIE/Mechanics/single_joint_displacement_controlled.prj b/Tests/Data/LIE/Mechanics/single_joint_displacement_controlled.prj
index 86685345ab2d45790221ae3637daef8c810b48eb..b8deb674a9051eb709b1a9116c92c86a83988bef 100644
--- a/Tests/Data/LIE/Mechanics/single_joint_displacement_controlled.prj
+++ b/Tests/Data/LIE/Mechanics/single_joint_displacement_controlled.prj
@@ -26,6 +26,10 @@
<cohesion>cohesion</cohesion>
<penalty_aperture_cutoff>1e-7</penalty_aperture_cutoff>
<tension_cutoff>1</tension_cutoff>
+ <nonlinear_solver>
+ <maximum_iterations>1000</maximum_iterations>
+ <error_tolerance>1e-12</error_tolerance>
+ </nonlinear_solver>
</fracture_model>
<fracture_properties>
<material_id>1</material_id>
@@ -43,9 +47,9 @@
<process ref="SD">
<nonlinear_solver>basic_newton</nonlinear_solver>
<convergence_criterion>
- <type>Residual</type>
+ <type>PerComponentDeltaX</type>
<norm_type>NORM2</norm_type>
- <abstol>1e-11</abstol>
+ <reltols>1e-15 1e-15 1e-15 1e-15 1e-15 1e-15</reltols>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
@@ -199,9 +203,7 @@
<lis>-i bicgstab -p jacobi -tol 1e-16 -maxiter 10000</lis>
<eigen>
<solver_type>SparseLU</solver_type>
- <precon_type>DIAGONAL</precon_type>
- <max_iteration_step>10000</max_iteration_step>
- <error_tolerance>1e-16</error_tolerance>
+ <scaling>1</scaling>
</eigen>
<petsc>
<prefix>sd</prefix>
diff --git a/Tests/MaterialLib/TestFractureModels.cpp b/Tests/MaterialLib/TestFractureModels.cpp
index 38690268c5c52b1d8d2cb94df04cab3403db19ed..7eaf274b2ce840e7dff20fd31c604f03d9d8d99f 100644
--- a/Tests/MaterialLib/TestFractureModels.cpp
+++ b/Tests/MaterialLib/TestFractureModels.cpp
@@ -21,8 +21,15 @@
using namespace MaterialLib::Fracture;
-static const double eps_sigma = 1e6*1e-5;
-static const double eps_C = 1e10*1e-5;
+// For known exact quantities like zero.
+constexpr double eps = std::numeric_limits<double>::epsilon();
+// For numeric in-order-of-sigma comparisons.
+constexpr double eps_sigma = 1e6 * eps;
+// For numeric in-order-of-C comparisons.
+constexpr double eps_C = 2e10 * eps;
+
+static NumLib::NewtonRaphsonSolverParameters const nonlinear_solver_parameters{
+ 1000, 1e-16};
TEST(MaterialLib_Fracture, LinearElasticIsotropic)
{
@@ -37,6 +44,7 @@ TEST(MaterialLib_Fracture, LinearElasticIsotropic)
tension_cutoff, mp};
std::unique_ptr<FractureModelBase<2>::MaterialStateVariables> state(
fractureModel.createMaterialStateVariables());
+ state->pushBackState();
Eigen::Vector2d const w_prev = Eigen::Vector2d::Zero();
Eigen::Vector2d const sigma0 = Eigen::Vector2d::Zero();
@@ -51,24 +59,23 @@ TEST(MaterialLib_Fracture, LinearElasticIsotropic)
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- 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);
-
+ EXPECT_NEAR(-1e4, sigma[0], eps_sigma);
+ EXPECT_NEAR(-1e6, sigma[1], eps_sigma);
+ EXPECT_NEAR(1e9, C(0, 0), eps);
+ EXPECT_NEAR(0, C(0, 1), eps);
+ EXPECT_NEAR(0, C(1, 0), eps);
+ EXPECT_NEAR(1e11, C(1, 1), eps);
w << -1e-5, 1e-5;
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- 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);
+ EXPECT_NEAR(0, sigma[0], eps);
+ EXPECT_NEAR(0, sigma[1], eps);
+ EXPECT_NEAR(0, C(0, 0), eps);
+ EXPECT_NEAR(0, C(0, 1), eps);
+ EXPECT_NEAR(0, C(1, 0), eps);
+ EXPECT_NEAR(0, C(1, 1), eps);
}
TEST(MaterialLib_Fracture, MohrCoulomb2D_elastic)
@@ -84,10 +91,12 @@ TEST(MaterialLib_Fracture, MohrCoulomb2D_elastic)
double const aperture0 = 1;
double const penalty_aperture_cutoff = aperture0;
bool const tension_cutoff = true;
- MohrCoulomb::MohrCoulomb<2> fractureModel{penalty_aperture_cutoff,
+ MohrCoulomb::MohrCoulomb<2> fractureModel{nonlinear_solver_parameters,
+ penalty_aperture_cutoff,
tension_cutoff, mp};
std::unique_ptr<FractureModelBase<2>::MaterialStateVariables> state(
fractureModel.createMaterialStateVariables());
+ state->pushBackState();
Eigen::Vector2d const w_prev = Eigen::Vector2d::Zero();
Eigen::Vector2d const sigma0 = Eigen::Vector2d::Zero();
@@ -102,23 +111,24 @@ TEST(MaterialLib_Fracture, MohrCoulomb2D_elastic)
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- 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);
+ EXPECT_NEAR(-1e4, sigma[0], eps_sigma);
+ EXPECT_NEAR(-1e6, sigma[1], eps_sigma);
+ EXPECT_NEAR(1e9, C(0, 0), eps);
+ EXPECT_NEAR(0, C(0, 1), eps);
+ EXPECT_NEAR(0, C(1, 0), eps);
+ EXPECT_NEAR(1e11, C(1, 1), eps);
w << -1e-5, 1e-5;
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- 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);
+ EXPECT_NEAR(0, sigma[0], eps);
+ EXPECT_NEAR(0, sigma[1], eps);
+ EXPECT_NEAR(0, C(0, 0), eps);
+ EXPECT_NEAR(0, C(0, 1), eps);
+ EXPECT_NEAR(0, C(1, 0), eps);
+ EXPECT_NEAR(0, C(1, 1), eps);
+ EXPECT_LE(state->getShearYieldFunctionValue(), 0);
}
TEST(MaterialLib_Fracture, MohrCoulomb2D_negative_t)
@@ -134,10 +144,12 @@ TEST(MaterialLib_Fracture, MohrCoulomb2D_negative_t)
double const aperture0 = 1e-5;
double const penalty_aperture_cutoff = aperture0;
bool const tension_cutoff = true;
- MohrCoulomb::MohrCoulomb<2> fractureModel{penalty_aperture_cutoff,
+ MohrCoulomb::MohrCoulomb<2> fractureModel{nonlinear_solver_parameters,
+ penalty_aperture_cutoff,
tension_cutoff, mp};
std::unique_ptr<FractureModelBase<2>::MaterialStateVariables> state(
fractureModel.createMaterialStateVariables());
+ state->pushBackState();
Eigen::Vector2d const w_prev = Eigen::Vector2d::Zero();
Eigen::Vector2d const sigma0(-3.46e6, -2e6);
@@ -152,13 +164,13 @@ TEST(MaterialLib_Fracture, MohrCoulomb2D_negative_t)
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- 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);
- ASSERT_NEAR(1.06608E+5, state->getShearYieldFunctionValue(), eps_sigma);
+ EXPECT_NEAR(-3575294.0369758257, sigma[0], eps_sigma);
+ EXPECT_NEAR(-2147026.5752851903, sigma[1], eps_sigma);
+ EXPECT_NEAR(1107234904.9501991, C(0, 0), eps_C);
+ EXPECT_NEAR(12655752875.023743, C(0, 1), eps_C);
+ EXPECT_NEAR(4132256921.1878347, C(1, 0), eps_C);
+ EXPECT_NEAR(47231912737.624504, C(1, 1), eps_C);
+ EXPECT_NEAR(0, state->getShearYieldFunctionValue(), eps);
}
@@ -175,10 +187,12 @@ TEST(MaterialLib_Fracture, MohrCoulomb2D_positive_t)
double const aperture0 = 1e-5;
double const penalty_aperture_cutoff = aperture0;
bool const tension_cutoff = true;
- MohrCoulomb::MohrCoulomb<2> fractureModel{penalty_aperture_cutoff,
+ MohrCoulomb::MohrCoulomb<2> fractureModel{nonlinear_solver_parameters,
+ penalty_aperture_cutoff,
tension_cutoff, mp};
std::unique_ptr<FractureModelBase<2>::MaterialStateVariables> state(
fractureModel.createMaterialStateVariables());
+ state->pushBackState();
Eigen::Vector2d const w_prev = Eigen::Vector2d::Zero();
Eigen::Vector2d const sigma0(0, -2e6);
@@ -193,13 +207,13 @@ TEST(MaterialLib_Fracture, MohrCoulomb2D_positive_t)
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- ASSERT_NEAR(253972.0, sigma[0], eps_sigma);
- ASSERT_NEAR(-2.94784e+06, 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);
- ASSERT_NEAR(446608.0, state->getShearYieldFunctionValue(), eps_sigma);
+ EXPECT_NEAR(3594117.0303599793, sigma[0], eps_sigma);
+ EXPECT_NEAR(-2217274.9429453835, sigma[1], eps_sigma);
+ EXPECT_NEAR(1107234904.9501991, C(0, 0), eps_C);
+ EXPECT_NEAR(-12655752875.023743, C(0, 1), eps_C);
+ EXPECT_NEAR(-4132256921.1878347, C(1, 0), eps_C);
+ EXPECT_NEAR(47231912737.624504, C(1, 1), eps_C);
+ EXPECT_NEAR(0, state->getShearYieldFunctionValue(), eps);
}
TEST(MaterialLib_Fracture, MohrCoulomb3D_negative_t1)
@@ -215,10 +229,12 @@ TEST(MaterialLib_Fracture, MohrCoulomb3D_negative_t1)
double const aperture0 = 1e-5;
double const penalty_aperture_cutoff = aperture0;
bool const tension_cutoff = true;
- MohrCoulomb::MohrCoulomb<3> fractureModel{penalty_aperture_cutoff,
+ MohrCoulomb::MohrCoulomb<3> fractureModel{nonlinear_solver_parameters,
+ penalty_aperture_cutoff,
tension_cutoff, mp};
std::unique_ptr<FractureModelBase<3>::MaterialStateVariables> state(
fractureModel.createMaterialStateVariables());
+ state->pushBackState();
Eigen::Vector3d const w_prev = Eigen::Vector3d::Zero();
Eigen::Vector3d const sigma0(-3.46e6, 0, -2e6);
@@ -233,19 +249,19 @@ TEST(MaterialLib_Fracture, MohrCoulomb3D_negative_t1)
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- ASSERT_NEAR(-3.50360e6, sigma[0], eps_sigma);
- ASSERT_NEAR(0.0, sigma[1], eps_sigma);
- ASSERT_NEAR(-2.16271e6, sigma[2], eps_sigma);
- ASSERT_NEAR(1.10723e+09, C(0,0), eps_C);
- ASSERT_NEAR(0.0, C(0,1), eps_C);
- ASSERT_NEAR(1.26558e+10, C(0,2), eps_C);
- ASSERT_NEAR(0.0, C(1,0), eps_C);
- ASSERT_NEAR(20.e9, C(1,1), eps_C);
- ASSERT_NEAR(0.0, C(1,2), eps_C);
- ASSERT_NEAR(4.13226e+09, C(2,0), eps_C);
- ASSERT_NEAR(0.0, C(2,1), eps_C);
- ASSERT_NEAR(4.72319e+10, C(2,2), eps_C);
- ASSERT_NEAR(1.06608E+5, state->getShearYieldFunctionValue(), eps_sigma);
+ EXPECT_NEAR(-3575294.0369758257, sigma[0], eps_sigma);
+ EXPECT_NEAR(0.0, sigma[1], eps);
+ EXPECT_NEAR(-2147026.5752851903, sigma[2], eps_sigma);
+ EXPECT_NEAR(1107234904.9501991, C(0, 0), eps_C);
+ EXPECT_NEAR(0.0, C(0, 1), eps);
+ EXPECT_NEAR(12655752875.023743, C(0, 2), eps_C);
+ EXPECT_NEAR(0.0, C(1, 0), eps);
+ EXPECT_NEAR(20.e9, C(1, 1), eps);
+ EXPECT_NEAR(0.0, C(1, 2), eps);
+ EXPECT_NEAR(4132256921.1878347, C(2, 0), eps_C);
+ EXPECT_NEAR(0.0, C(2, 1), eps);
+ EXPECT_NEAR(47231912737.624504, C(2, 2), eps_C);
+ EXPECT_NEAR(0, state->getShearYieldFunctionValue(), eps);
}
TEST(MaterialLib_Fracture, MohrCoulomb3D_positive_t1)
@@ -261,10 +277,12 @@ TEST(MaterialLib_Fracture, MohrCoulomb3D_positive_t1)
double const aperture0 = 1e-5;
double const penalty_aperture_cutoff = aperture0;
bool const tension_cutoff = true;
- MohrCoulomb::MohrCoulomb<3> fractureModel{penalty_aperture_cutoff,
+ MohrCoulomb::MohrCoulomb<3> fractureModel{nonlinear_solver_parameters,
+ penalty_aperture_cutoff,
tension_cutoff, mp};
std::unique_ptr<FractureModelBase<3>::MaterialStateVariables> state(
fractureModel.createMaterialStateVariables());
+ state->pushBackState();
Eigen::Vector3d const w_prev = Eigen::Vector3d::Zero();
Eigen::Vector3d const sigma0(0, 0, -2e6);
@@ -279,19 +297,19 @@ TEST(MaterialLib_Fracture, MohrCoulomb3D_positive_t1)
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- ASSERT_NEAR(253972.0, sigma[0], eps_sigma);
- ASSERT_NEAR(0.0, sigma[1], eps_sigma);
- ASSERT_NEAR(-2.94784e+06, sigma[2], eps_sigma);
- ASSERT_NEAR(1.10723e+09, C(0,0), eps_C);
- ASSERT_NEAR(0.0, C(0,1), eps_C);
- ASSERT_NEAR(-1.26558e+10, C(0,2), eps_C);
- ASSERT_NEAR(0.0, C(1,0), eps_C);
- ASSERT_NEAR(20.e9, C(1,1), eps_C);
- ASSERT_NEAR(0.0, C(1,2), eps_C);
- ASSERT_NEAR(-4.13226e+09, C(2,0), eps_C);
- ASSERT_NEAR(0.0, C(2,1), eps_C);
- ASSERT_NEAR(4.72319e+10, C(2,2), eps_C);
- ASSERT_NEAR(446608.0, state->getShearYieldFunctionValue(), eps_sigma);
+ EXPECT_NEAR(3594117.0303599793, sigma[0], eps_sigma);
+ EXPECT_NEAR(0.0, sigma[1], eps);
+ EXPECT_NEAR(-2217274.9429453835, sigma[2], eps_sigma);
+ EXPECT_NEAR(1107234904.9501991, C(0, 0), eps_C);
+ EXPECT_NEAR(0.0, C(0, 1), eps);
+ EXPECT_NEAR(-12655752875.023743, C(0, 2), eps_C);
+ EXPECT_NEAR(0.0, C(1, 0), eps);
+ EXPECT_NEAR(20.e9, C(1, 1), eps);
+ EXPECT_NEAR(0.0, C(1, 2), eps);
+ EXPECT_NEAR(-4132256921.1878347, C(2, 0), eps_C);
+ EXPECT_NEAR(0.0, C(2, 1), eps);
+ EXPECT_NEAR(47231912737.624504, C(2, 2), eps_C);
+ EXPECT_NEAR(0, state->getShearYieldFunctionValue(), eps);
}
@@ -308,10 +326,12 @@ TEST(MaterialLib_Fracture, MohrCoulomb3D_positive_t2)
double const aperture0 = 1e-5;
double const penalty_aperture_cutoff = aperture0;
bool const tension_cutoff = true;
- MohrCoulomb::MohrCoulomb<3> fractureModel{penalty_aperture_cutoff,
+ MohrCoulomb::MohrCoulomb<3> fractureModel{nonlinear_solver_parameters,
+ penalty_aperture_cutoff,
tension_cutoff, mp};
std::unique_ptr<FractureModelBase<3>::MaterialStateVariables> state(
fractureModel.createMaterialStateVariables());
+ state->pushBackState();
Eigen::Vector3d const w_prev = Eigen::Vector3d::Zero();
Eigen::Vector3d const sigma0(0, 0, -2e6);
@@ -326,19 +346,19 @@ TEST(MaterialLib_Fracture, MohrCoulomb3D_positive_t2)
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- ASSERT_NEAR(0.0, sigma[0], eps_sigma);
- ASSERT_NEAR(253972.0, sigma[1], eps_sigma);
- ASSERT_NEAR(-2.94784e+06, sigma[2], eps_sigma);
- ASSERT_NEAR(20.e9, C(0,0), eps_C);
- ASSERT_NEAR(0.0, C(0,1), eps_C);
- ASSERT_NEAR(0.0, C(0,2), eps_C);
- ASSERT_NEAR(0.0, C(1,0), eps_C);
- ASSERT_NEAR(1.10723e+09, C(1,1), eps_C);
- ASSERT_NEAR(-1.26558e+10, C(1,2), eps_C);
- ASSERT_NEAR(0.0, C(2,0), eps_C);
- ASSERT_NEAR(-4.13226e+09, C(2,1), eps_C);
- ASSERT_NEAR(4.72319e+10, C(2,2), eps_C);
- ASSERT_NEAR(446608.0, state->getShearYieldFunctionValue(), eps_sigma);
+ EXPECT_NEAR(0.0, sigma[0], eps);
+ EXPECT_NEAR(3594117.0303599793, sigma[1], eps_sigma);
+ EXPECT_NEAR(-2217274.9429453835, sigma[2], eps_sigma);
+ EXPECT_NEAR(20.e9, C(0, 0), eps);
+ EXPECT_NEAR(0.0, C(0, 1), eps);
+ EXPECT_NEAR(0.0, C(0, 2), eps);
+ EXPECT_NEAR(0.0, C(1, 0), eps);
+ EXPECT_NEAR(1107234904.9501991, C(1, 1), eps_C);
+ EXPECT_NEAR(-12655752875.023743, C(1, 2), eps_C);
+ EXPECT_NEAR(0.0, C(2, 0), eps);
+ EXPECT_NEAR(-4132256921.1878347, C(2, 1), eps_C);
+ EXPECT_NEAR(47231912737.624504, C(2, 2), eps_C);
+ EXPECT_NEAR(0, state->getShearYieldFunctionValue(), eps);
}
@@ -355,10 +375,11 @@ TEST(MaterialLib_Fracture, MohrCoulomb3D_negative_t1t2)
double const aperture0 = 1e-5;
double const penalty_aperture_cutoff = aperture0;
bool const tension_cutoff = true;
- MohrCoulomb::MohrCoulomb<3> fractureModel{penalty_aperture_cutoff,
- tension_cutoff, mp};
+ MohrCoulomb::MohrCoulomb<3> fractureModel{
+ {1000, 1e-9}, penalty_aperture_cutoff, tension_cutoff, mp};
std::unique_ptr<FractureModelBase<3>::MaterialStateVariables> state(
fractureModel.createMaterialStateVariables());
+ state->pushBackState();
Eigen::Vector3d const w_prev = Eigen::Vector3d::Zero();
Eigen::Vector3d const sigma0(-3.46e6 / std::sqrt(2), -3.46e6 / std::sqrt(2),
@@ -374,19 +395,20 @@ TEST(MaterialLib_Fracture, MohrCoulomb3D_negative_t1t2)
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- ASSERT_NEAR(-3.50360e6/std::sqrt(2), sigma[0], eps_sigma);
- ASSERT_NEAR(-3.50360e6/std::sqrt(2), sigma[1], eps_sigma);
- ASSERT_NEAR(-2.16271e6, sigma[2], eps_sigma);
- ASSERT_NEAR(1.05536e+10, C(0,0), eps_C);
- ASSERT_NEAR(-9.44638e+09, C(0,1), eps_C);
- ASSERT_NEAR(8.94897e+09, C(0,2), eps_C);
- ASSERT_NEAR(-9.44638e+09, C(1,0), eps_C);
- ASSERT_NEAR(1.05536e+10, C(1,1), eps_C);
- ASSERT_NEAR(8.94897e+09, C(1,2), eps_C);
- ASSERT_NEAR(2.92195e+09, C(2,0), eps_C);
- ASSERT_NEAR(2.92195e+09, C(2,1), eps_C);
- ASSERT_NEAR(4.72319e+10, C(2,2), eps_C);
- ASSERT_NEAR(1.06608E+5, state->getShearYieldFunctionValue(), eps_sigma);
+ EXPECT_NEAR(-3575294.0369758265 / std::sqrt(2), sigma[0], eps_sigma);
+ EXPECT_NEAR(-3575294.0369758265 / std::sqrt(2), sigma[1], eps_sigma);
+ EXPECT_NEAR(-2147026.5752851903, sigma[2], eps_sigma);
+ EXPECT_NEAR(10553617452.475098, C(0, 0), eps_C);
+ EXPECT_NEAR(-9446382547.5249023, C(0, 1), eps_C);
+ EXPECT_NEAR(8948968678.9504356, C(0, 2), eps_C);
+ EXPECT_NEAR(-9446382547.5249023, C(1, 0), eps_C);
+ EXPECT_NEAR(10553617452.475098, C(1, 1), eps_C);
+ EXPECT_NEAR(8948968678.9504356, C(1, 2), eps_C);
+ EXPECT_NEAR(2921946890.5769634, C(2, 0), eps_C);
+ EXPECT_NEAR(2921946890.5769634, C(2, 1), eps_C);
+ EXPECT_NEAR(47231912737.624504, C(2, 2), eps_C);
+ EXPECT_NEAR(0, state->getShearYieldFunctionValue(),
+ 1e-9); // same as newton tolerance
}
@@ -403,10 +425,12 @@ TEST(MaterialLib_Fracture, MohrCoulomb3D_negative_t1_positive_t2)
double const aperture0 = 1e-5;
double const penalty_aperture_cutoff = aperture0;
bool const tension_cutoff = true;
- MohrCoulomb::MohrCoulomb<3> fractureModel{penalty_aperture_cutoff,
- tension_cutoff, mp};
+
+ MohrCoulomb::MohrCoulomb<3> fractureModel{
+ {1000, 1e-9}, penalty_aperture_cutoff, tension_cutoff, mp};
std::unique_ptr<FractureModelBase<3>::MaterialStateVariables> state(
fractureModel.createMaterialStateVariables());
+ state->pushBackState();
Eigen::Vector3d const w_prev = Eigen::Vector3d::Zero();
Eigen::Vector3d const sigma0(-3.46e6 / std::sqrt(2), 3.46e6 / std::sqrt(2),
@@ -422,18 +446,19 @@ TEST(MaterialLib_Fracture, MohrCoulomb3D_negative_t1_positive_t2)
fractureModel.computeConstitutiveRelation(0, x, aperture0, sigma0, w_prev,
w, sigma_prev, sigma, C, *state);
- ASSERT_NEAR(-3.50360e6/std::sqrt(2), sigma[0], eps_sigma);
- ASSERT_NEAR(3.50360e6/std::sqrt(2), sigma[1], eps_sigma);
- ASSERT_NEAR(-2.16271e6, sigma[2], eps_sigma);
- ASSERT_NEAR(1.05536e+10, C(0,0), eps_C);
- ASSERT_NEAR(9.44638e+09, C(0,1), eps_C);
- ASSERT_NEAR(8.94897e+09, C(0,2), eps_C);
- ASSERT_NEAR(9.44638e+09, C(1,0), eps_C);
- ASSERT_NEAR(1.05536e+10, C(1,1), eps_C);
- ASSERT_NEAR(-8.94897e+09, C(1,2), eps_C);
- ASSERT_NEAR(2.92195e+09, C(2,0), eps_C);
- ASSERT_NEAR(-2.92195e+09, C(2,1), eps_C);
- ASSERT_NEAR(4.72319e+10, C(2,2), eps_C);
- ASSERT_NEAR(1.06608E+5, state->getShearYieldFunctionValue(), eps_sigma);
+ EXPECT_NEAR(-3575294.0369758265 / std::sqrt(2), sigma[0], eps_sigma);
+ EXPECT_NEAR(3575294.0369758265 / std::sqrt(2), sigma[1], eps_sigma);
+ EXPECT_NEAR(-2147026.5752851903, sigma[2], eps_sigma);
+ EXPECT_NEAR(10553617452.475098, C(0, 0), eps_C);
+ EXPECT_NEAR(9446382547.5249023, C(0, 1), eps_C);
+ EXPECT_NEAR(8948968678.9504356, C(0, 2), eps_C);
+ EXPECT_NEAR(9446382547.5249023, C(1, 0), eps_C);
+ EXPECT_NEAR(10553617452.475098, C(1, 1), eps_C);
+ EXPECT_NEAR(-8948968678.9504356, C(1, 2), eps_C);
+ EXPECT_NEAR(2921946890.5769634, C(2, 0), eps_C);
+ EXPECT_NEAR(-2921946890.5769634, C(2, 1), eps_C);
+ EXPECT_NEAR(47231912737.624504, C(2, 2), eps_C);
+ EXPECT_NEAR(0, state->getShearYieldFunctionValue(),
+ 1e-9); // same as newton tolerance
}