diff --git a/MaterialLib/SolidModels/MFront/CMakeLists.txt b/MaterialLib/SolidModels/MFront/CMakeLists.txt
index b1b5fa5c1ae5c8d91d1acf890d68bfc1eed07894..2a9678ef3657d0f82f27ee1b95afd5cae5301e51 100644
--- a/MaterialLib/SolidModels/MFront/CMakeLists.txt
+++ b/MaterialLib/SolidModels/MFront/CMakeLists.txt
@@ -20,7 +20,10 @@ mfront_behaviours_check_library(OgsMFrontBehaviour
StandardElasticityBrick
StandardElasticityBrickOrtho
Lubby2
- Lubby2mod)
+ Lubby2mod
+ MohrCoulombAbboSloanAniso
+ MohrCoulombAbboSloanUBI
+ MohrCoulombAbboSloanUBIOrtho)
target_include_directories(MaterialLib_SolidModels_MFront
PUBLIC ThirdParty/MGIS/include)
diff --git a/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloan.mfront b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloan.mfront
index e42b007fbec6a2be90c143a48f62f319feaabae1..36e6bb528063f21f9b6ab45618c98e1f36285a8f 100644
--- a/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloan.mfront
+++ b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloan.mfront
@@ -85,7 +85,7 @@ a.setEntryName("TensionCutOffParameter");
// Compute initial elastic strain
const auto S = invert(D);
- eel = S*sig;
+ eel = S * sig;
// elastic prediction
const auto sig_el = computeElasticPrediction();
@@ -116,7 +116,7 @@ a.setEntryName("TensionCutOffParameter");
}
const auto sMC =
I1_el / 3 * sin_phi + sqrt(J2_el * K * K + a * a * sin_phi * sin_phi);
- F = sMC - c * cos_phi > 0.;
+ F = sMC - c * cos_phi > 0.;
np = Stensor(real(0));
}
@@ -190,7 +190,8 @@ a.setEntryName("TensionCutOffParameter");
const auto dev_s_squared = computeJ3Derivative(
sig); // replaces dev_s_squared = deviator(square(s));
const auto dG_dI1 = sin_psi / 3.;
- const auto root = max(sqrt(J2 * KG * KG + a * a * tan(phi) * tan(phi) * cos(psi) * cos(psi)),
+ const auto root = max(sqrt(J2 * KG * KG + a * a * tan(phi) * tan(phi) *
+ cos(psi) * cos(psi)),
local_zero_tolerance);
const auto dG_dJ2 = KG / (2. * root) * (KG - tan_3_lode * dKG_dlode);
const auto dG_dJ3 = J2 * KG * tan_3_lode / (3. * J3 * root) * dKG_dlode;
@@ -202,10 +203,9 @@ a.setEntryName("TensionCutOffParameter");
// yield function
const auto rootF = max(sqrt(J2 * K * K + a * a * sin_phi * sin_phi),
local_zero_tolerance);
- const auto Fy1 = I1 * sin_phi / 3 + rootF;
- const auto Fy = Fy1 - c * cos_phi;
+ const auto Fy = I1 * sin_phi / 3 + rootF - c * cos_phi;
// return if overshoot of yield surface
- if (Fy1 - c * cos_phi > 1e-4 * D(0, 0))
+ if (Fy > 1e-4 * D(0, 0))
{
return false;
}
@@ -217,7 +217,7 @@ a.setEntryName("TensionCutOffParameter");
const auto nF = eval(dF_dI1 * id + dF_dJ2 * s + dF_dJ3 * dev_s_squared);
// building dfeel_ddeel
- const auto Pdev = id4 - (id ^ id) / 3;
+ const auto Pdev = Stensor4::K();
const auto dG_dlode = KG * J2 / (root)*dKG_dlode;
const auto dG_ddlode =
diff --git a/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanAniso.mfront b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanAniso.mfront
new file mode 100644
index 0000000000000000000000000000000000000000..b4390b716dfb967bb3e3a7a8e1df48b21f2793c6
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanAniso.mfront
@@ -0,0 +1,285 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+@DSL Implicit;
+@Behaviour MohrCoulombAbboSloanAniso;
+@Author Thomas Nagel;
+@Date 21 / 03 / 2020;
+
+@Algorithm NewtonRaphson;
+@MaximumNumberOfIterations 200;
+
+@Description
+{
+Non-associated Mohr Coulomb model with singularity smoothing according to
+Abbo-Sloan.Anisotropy is introduced by stress-scaling according to
+//
+Malcom, M.A.M.(2018).Analysis of underground excavations in argillaceous hard
+soils - weak rocks. Technical University of Catalonia.
+//
+The plasticity is transversely isotropic; cn always scales the(1, 0, 0)
+direction. Rotations performed by the FE Programme.
+};
+
+@OrthotropicBehaviour<Pipe>;
+@Brick StandardElasticity;
+
+@Theta 1.0; // time integration scheme
+@Epsilon 1e-14; // tolerance of local stress integration algorithm
+@ModellingHypotheses{".+"};
+
+@RequireStiffnessTensor<UnAltered>;
+
+@StateVariable real lam;
+lam.setGlossaryName("EquivalentPlasticStrain");
+
+@Parameter pi = 3.14159265359;
+@Parameter local_zero_tolerance = 1.e-14;
+
+// Note: YoungModulus and PoissonRatio defined as parameters
+// Note: Glossary names are already given; entry names are newly defined
+
+@MaterialProperty stress c;
+c.setEntryName("Cohesion");
+@MaterialProperty real phi;
+phi.setEntryName("FrictionAngle");
+@MaterialProperty real psi;
+psi.setEntryName("DilatancyAngle");
+@MaterialProperty real lodeT;
+lodeT.setEntryName("TransitionAngle");
+@MaterialProperty stress a;
+a.setEntryName("TensionCutOffParameter");
+@MaterialProperty real cn;
+cn.setEntryName("NormalFactor");
+@MaterialProperty real cs;
+cs.setEntryName("ShearFactor");
+
+@LocalVariable Stensor np;
+
+@LocalVariable bool F; // if true, plastic loading
+@LocalVariable real sin_psi;
+@LocalVariable real sin_phi;
+@LocalVariable real cos_phi;
+@LocalVariable real cos_lodeT;
+@LocalVariable real sin_lodeT;
+@LocalVariable real tan_lodeT;
+@LocalVariable real cos_3_lodeT;
+@LocalVariable real sin_3_lodeT;
+@LocalVariable real tan_3_lodeT;
+@LocalVariable real multi;
+
+@InitLocalVariables
+{
+ constexpr auto sqrt3 = Cste<real>::sqrt3;
+ constexpr auto isqrt3 = Cste<real>::isqrt3;
+ // conversion to rad
+ phi *= pi / 180.;
+ psi *= pi / 180.;
+ lodeT *= pi / 180.;
+ sin_psi = sin(psi);
+ cos_phi = cos(phi);
+ sin_phi = sin(phi);
+ sin_lodeT = sin(lodeT);
+ cos_lodeT = cos(lodeT);
+ tan_lodeT = tan(lodeT);
+ cos_3_lodeT = cos(3. * lodeT);
+ sin_3_lodeT = sin(3. * lodeT);
+ tan_3_lodeT = tan(3. * lodeT);
+
+ auto Hill = Stensor4::Id();
+ const int step = sig.size() + 1;
+ Hill[0] = cn;
+ Hill[step] = Hill[step * 2] = 1. / cn;
+ Hill[step * 3] = Hill[step * 4] = cs;
+
+ // Compute initial elastic strain
+ const auto S = invert(D);
+ eel = S * sig;
+
+ // elastic prediction
+ const auto sig_el = computeElasticPrediction();
+ const auto sig_el_scaled = Hill * sig_el;
+ const auto s_el_scaled = deviator(sig_el_scaled);
+ const auto I1_el = trace(sig_el_scaled);
+ const auto J2_el =
+ max((s_el_scaled | s_el_scaled) / 2., local_zero_tolerance);
+ const auto J3_el = det(s_el_scaled);
+ const auto arg = min(max(-3. * sqrt3 * J3_el / (2. * J2_el * sqrt(J2_el)),
+ -1. + local_zero_tolerance),
+ 1. - local_zero_tolerance);
+ const auto lode_el = 1. / 3. * asin(arg);
+ auto K = 0.0;
+ if (fabs(lode_el) < lodeT)
+ {
+ K = cos(lode_el) - isqrt3 * sin_phi * sin(lode_el);
+ }
+ else
+ {
+ const auto sign = min(
+ max(lode_el / max(abs(lode_el), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_phi);
+ const auto B = 1 / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_phi * cos_lodeT);
+ K = A - B * arg;
+ }
+ const auto sMC =
+ I1_el / 3 * sin_phi + sqrt(J2_el * K * K + a * a * sin_phi * sin_phi);
+ F = sMC - c * cos_phi > 0.;
+ np = Stensor(real(0));
+}
+
+@Integrator
+{
+ constexpr auto sqrt3 = Cste<real>::sqrt3;
+ constexpr auto isqrt3 = Cste<real>::isqrt3;
+ constexpr auto id = Stensor::Id();
+ auto Hill = Stensor4::Id();
+ const int step = sig.size() + 1;
+ Hill[0] = cn;
+ Hill[step] = Hill[step * 2] = 1. / cn;
+ Hill[step * 3] = Hill[step * 4] = cs;
+ if (F)
+ {
+ const auto sig_scaled = Hill * sig;
+ const auto s = deviator(sig_scaled);
+ const auto I1 = trace(sig_scaled);
+ const auto J2 = max((s | s) / 2., local_zero_tolerance);
+ const auto J3 = real(det(s) < 0. ? min(det(s), -local_zero_tolerance)
+ : max(det(s), local_zero_tolerance));
+ const auto arg = min(max(-3. * sqrt3 * J3 / (2. * J2 * sqrt(J2)),
+ -1. + local_zero_tolerance),
+ 1. - local_zero_tolerance);
+ const auto lode = 1. / 3. * asin(arg);
+ const auto cos_lode = cos(lode);
+ const auto sin_lode = sin(lode);
+ const auto cos_3_lode = cos(3. * lode);
+ const auto sin_3_lode = arg;
+ const auto tan_3_lode = tan(3. * lode);
+ auto K = 0.;
+ auto dK_dlode = 1.;
+ if (fabs(lode) < lodeT)
+ {
+ K = cos_lode - isqrt3 * sin_phi * sin_lode;
+ dK_dlode = -sin_lode - isqrt3 * sin_phi * cos_lode;
+ }
+ else
+ {
+ const auto sign =
+ min(max(lode / max(abs(lode), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_phi);
+ const auto B = 1. / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_phi * cos_lodeT);
+ K = A - B * sin_3_lode;
+ dK_dlode = -3. * B * cos_3_lode;
+ }
+ auto KG = 0.0; // move into a function to avoid code duplication
+ auto dKG_dlode = 1.;
+ auto dKG_ddlode = 1.;
+ if (fabs(lode) < lodeT)
+ {
+ KG = cos_lode - isqrt3 * sin_psi * sin_lode;
+ dKG_dlode = -sin_lode - isqrt3 * sin_psi * cos_lode;
+ dKG_ddlode = -cos_lode + isqrt3 * sin_psi * sin_lode;
+ }
+ else
+ {
+ const auto sign =
+ min(max(lode / max(fabs(lode), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_psi);
+ const auto B = 1. / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_psi * cos_lodeT);
+ KG = A - B * sin_3_lode;
+ dKG_dlode = -3. * B * cos_3_lode;
+ dKG_ddlode = 9. * B * sin_3_lode;
+ }
+
+ // flow direction
+ const auto dev_s_squared = computeJ3Derivative(
+ sig_scaled); // replaces dev_s_squared = deviator(square(s));
+ const auto dG_dI1 = sin_psi / 3.;
+ const auto root = max(sqrt(J2 * KG * KG + a * a * tan(phi) * tan(phi) *
+ cos(psi) * cos(psi)),
+ local_zero_tolerance);
+ const auto dG_dJ2 = KG / (2. * root) * (KG - tan_3_lode * dKG_dlode);
+ const auto dG_dJ3 = J2 * KG * tan_3_lode / (3. * J3 * root) * dKG_dlode;
+ const auto n =
+ eval(Hill * (dG_dI1 * id + dG_dJ2 * s + dG_dJ3 * dev_s_squared));
+ if (this->iter > 30 && abs(n | np) < sqrt(n | n) * sqrt(np | np) * 0.99)
+ {
+ return false;
+ }
+ // yield function
+ const auto rootF = max(sqrt(J2 * K * K + a * a * sin_phi * sin_phi),
+ local_zero_tolerance);
+ const auto Fy = I1 * sin_phi / 3 + rootF - c * cos_phi;
+ // return if overshoot of yield surface
+ if (Fy > 1e-4 * D(0, 0))
+ {
+ return false;
+ }
+
+ // yield function derivative for Jacobian
+ const auto dF_dI1 = sin_phi / 3.;
+ const auto dF_dJ2 = K / (2. * rootF) * (K - tan_3_lode * dK_dlode);
+ const auto dF_dJ3 = J2 * K * tan_3_lode / (3. * J3 * rootF) * dK_dlode;
+ const auto nF =
+ eval(Hill * (dF_dI1 * id + dF_dJ2 * s + dF_dJ3 * dev_s_squared));
+
+ // building dfeel_ddeel
+ const auto Pdev = Stensor4::K();
+
+ const auto dG_dlode = KG * J2 / (root)*dKG_dlode;
+ const auto dG_ddlode =
+ J2 / root *
+ (dKG_dlode * dKG_dlode * (1. - J2 * KG * KG / (root * root)) +
+ KG * dKG_ddlode);
+ const auto dG_ddlodeJ2 =
+ KG / root * dKG_dlode * (1. - J2 * KG * KG / (2 * root * root));
+ const auto dG_ddJ2 =
+ -KG * KG * KG * KG / (4. * root * root * root) +
+ dG_dlode * tan_3_lode / (2 * J2 * J2) -
+ tan_3_lode / (2 * J2) *
+ (2 * dG_ddlodeJ2 - tan_3_lode / (2 * J2) * dG_ddlode -
+ 3 / (2 * J2 * cos_3_lode * cos_3_lode) * dG_dlode);
+ const auto dG_ddJ3 =
+ -tan_3_lode / (3 * J3 * J3) * dG_dlode +
+ tan_3_lode / (3 * J3) *
+ (dG_ddlode * tan_3_lode / (3 * J3) +
+ dG_dlode * 1. / (J3 * cos_3_lode * cos_3_lode));
+ const auto dG_ddJ2J3 =
+ dG_ddlodeJ2 * tan_3_lode / (3 * J3) -
+ tan_3_lode / (2 * J2) *
+ (dG_ddlode * tan_3_lode / (3 * J3) +
+ dG_dlode * 1. / (J3 * cos_3_lode * cos_3_lode));
+
+ // elasticity
+ feel += dlam * n;
+ dfeel_ddeel +=
+ theta * dlam * Hill *
+ (dG_dJ2 * Pdev + dG_dJ3 * computeJ3SecondDerivative(sig_scaled) +
+ dG_ddJ2 * (s ^ s) + dG_ddJ3 * (dev_s_squared ^ dev_s_squared) +
+ dG_ddJ2J3 * ((dev_s_squared ^ s) + (s ^ dev_s_squared))) *
+ Hill * D;
+ dfeel_ddlam = n;
+ // plasticity
+ flam = Fy / D(0, 0);
+ dflam_ddlam = strain(0);
+ dflam_ddeel = theta * (nF | D) / D(0, 0);
+ np = n;
+ }
+}
diff --git a/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanOrtho.mfront b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanOrtho.mfront
index c5f1aa7aefdcf2e5dc37d6fd3cb560337cda1a2c..c736c13fe141c35b2cb255403e61b2b7945b0302 100644
--- a/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanOrtho.mfront
+++ b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanOrtho.mfront
@@ -10,10 +10,11 @@
@DSL Implicit;
@Behaviour MohrCoulombAbboSloanOrtho;
@Author Thomas Nagel;
-@Date 05/02/2019;
+@Date 05 / 02 / 2019;
@Algorithm NewtonRaphson;
-//@Algorithm LevenbergMarquardt;
+@MaximumNumberOfIterations 200;
+// @Algorithm LevenbergMarquardt;
//@CompareToNumericalJacobian true;
//@NumericallyComputedJacobianBlocks {dfeel_ddeel};
// remove the above blocks once an analytical one is provided.
@@ -23,7 +24,7 @@
@OrthotropicBehaviour<Pipe>;
@Brick StandardElasticity;
-@Theta 1.0; // time integration scheme
+@Theta 1.0; // time integration scheme
@Epsilon 1e-14; // tolerance of local stress integration algorithm
@ModellingHypotheses{".+"};
@@ -83,7 +84,7 @@ a.setEntryName("TensionCutOffParameter");
// Compute initial elastic strain
const auto S = invert(D);
- eel = S*sig;
+ eel = S * sig;
// elastic prediction
const auto sig_el = computeElasticPrediction();
@@ -114,8 +115,7 @@ a.setEntryName("TensionCutOffParameter");
}
const auto sMC =
I1_el / 3 * sin_phi + sqrt(J2_el * K * K + a * a * sin_phi * sin_phi);
- //(int)(d * 1000.0)/1000.0
- F = (sMC - c * cos_phi) / D(0,0) > 0.;
+ F = sMC - c * cos_phi > 0.;
np = Stensor(real(0));
}
@@ -189,22 +189,22 @@ a.setEntryName("TensionCutOffParameter");
const auto dev_s_squared = computeJ3Derivative(
sig); // replaces dev_s_squared = deviator(square(s));
const auto dG_dI1 = sin_psi / 3.;
- const auto root = max(sqrt(J2 * KG * KG + a * a * sin_psi * sin_psi),
+ const auto root = max(sqrt(J2 * KG * KG + a * a * tan(phi) * tan(phi) *
+ cos(psi) * cos(psi)),
local_zero_tolerance);
const auto dG_dJ2 = KG / (2. * root) * (KG - tan_3_lode * dKG_dlode);
const auto dG_dJ3 = J2 * KG * tan_3_lode / (3. * J3 * root) * dKG_dlode;
const auto n = eval(dG_dI1 * id + dG_dJ2 * s + dG_dJ3 * dev_s_squared);
- if (this->iter > 30 && abs(n | np) < sqrt(n | n) * sqrt(np | np) * 0.99)
+ if (this->iter > 30 && abs(n | np) < sqrt(n | n) * sqrt(np | np) * 0.99)
{
return false;
}
- //yield function
+ // yield function
const auto rootF = max(sqrt(J2 * K * K + a * a * sin_phi * sin_phi),
local_zero_tolerance);
- const auto Fy1 = I1 * sin_phi / 3 + rootF;
- const auto Fy = Fy1 - c * cos_phi;
+ const auto Fy = I1 * sin_phi / 3 + rootF - c * cos_phi;
// return if overshoot of yield surface
- if (Fy1 - c * cos_phi > 1e-4 * D(0, 0))
+ if (Fy > 1e-4 * D(0, 0))
{
return false;
}
@@ -215,7 +215,7 @@ a.setEntryName("TensionCutOffParameter");
const auto nF = eval(dF_dI1 * id + dF_dJ2 * s + dF_dJ3 * dev_s_squared);
// building dfeel_ddeel
- const auto Pdev = id4 - (id ^ id) / 3;
+ const auto Pdev = Stensor4::K();
const auto dG_dlode = KG * J2 / (root)*dKG_dlode;
const auto dG_ddlode =
@@ -254,6 +254,6 @@ a.setEntryName("TensionCutOffParameter");
flam = Fy / D(0, 0);
dflam_ddlam = strain(0);
dflam_ddeel = theta * (nF | D) / D(0, 0);
- np = n;
+ np = n;
}
}
diff --git a/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanUBI.mfront b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanUBI.mfront
new file mode 100644
index 0000000000000000000000000000000000000000..b9582a5716fb2478958dd3b8cd32caf6d3f72341
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanUBI.mfront
@@ -0,0 +1,440 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+@DSL Implicit;
+@Behaviour MohrCoulombAbboSloanUBI;
+@Author Thomas Nagel;
+@Date 21 / 03 / 2020;
+@Description{
+Ubiquitous joint model where the matrix fails according to a Mohr-Coulomb
+criterion(with tension-cutoff) and the embedded weakness planes follow
+a Coulomb behaviour(so far without tension cut-off). Both models are
+non-associated.
+The fracture normal is the x-direction(1, 0, 0).
+Any rotations to be done by the FE programme.
+};
+
+@Algorithm NewtonRaphson;
+@MaximumNumberOfIterations 200;
+
+@Brick StandardElasticity;
+
+@Epsilon 1.e-14;
+@Theta 1.0;
+@Parameter local_zero_tolerance = 1.e-14;
+@Parameter pi = 3.14159265359;
+
+@ModellingHypotheses{".+"};
+@RequireStiffnessTensor<UnAltered>;
+
+// weakness plane cohesion
+@MaterialProperty real c_wp;
+c_wp.setEntryName("PlaneCohesion");
+// weakness plane friction angle
+@MaterialProperty real phi_wp;
+phi_wp.setEntryName("PlaneFrictionAngle");
+// weakness plane dilatancy angle
+@MaterialProperty real psi_wp;
+psi_wp.setEntryName("PlaneDilatancyAngle");
+
+@MaterialProperty stress c;
+c.setEntryName("Cohesion");
+@MaterialProperty real phi;
+phi.setEntryName("FrictionAngle");
+@MaterialProperty real psi;
+psi.setEntryName("DilatancyAngle");
+@MaterialProperty real lodeT;
+lodeT.setEntryName("TransitionAngle");
+@MaterialProperty stress a;
+a.setEntryName("TensionCutOffParameter");
+
+@LocalVariable Stensor np;
+@LocalVariable Stensor npwp;
+
+// Not array because of OGS output
+@StateVariable strain lamWP;
+lamWP.setEntryName("EquivalentPlasticStrainWP");
+@StateVariable strain lam;
+lam.setEntryName("EquivalentPlasticStrainMatrix");
+
+@LocalVariable bool F[2];
+@LocalVariable real tan_phi_wp;
+@LocalVariable real tan_psi_wp;
+@LocalVariable real sin_psi;
+@LocalVariable real sin_phi;
+@LocalVariable real cos_phi;
+@LocalVariable real cos_lodeT;
+@LocalVariable real sin_lodeT;
+@LocalVariable real tan_lodeT;
+@LocalVariable real cos_3_lodeT;
+@LocalVariable real sin_3_lodeT;
+@LocalVariable real tan_3_lodeT;
+
+@InitLocalVariables
+{
+ // tan_phi_wp after conversion to rad
+ tan_phi_wp = tan(phi_wp * pi / 180.);
+ // tan_psi_wp after conversion to rad
+ tan_psi_wp = tan(psi_wp * pi / 180.);
+
+ constexpr auto sqrt3 = Cste<real>::sqrt3;
+ constexpr auto isqrt3 = Cste<real>::isqrt3;
+ // conversion to rad
+ phi *= pi / 180.;
+ psi *= pi / 180.;
+ lodeT *= pi / 180.;
+ sin_psi = sin(psi);
+ cos_phi = cos(phi);
+ sin_phi = sin(phi);
+ sin_lodeT = sin(lodeT);
+ cos_lodeT = cos(lodeT);
+ tan_lodeT = tan(lodeT);
+ cos_3_lodeT = cos(3. * lodeT);
+ sin_3_lodeT = sin(3. * lodeT);
+ tan_3_lodeT = tan(3. * lodeT);
+
+ // Compute initial elastic strain
+ const auto S = invert(D);
+ eel = S * sig;
+
+ const auto sig_el = computeElasticPrediction();
+
+ const auto s_el = deviator(sig_el);
+ const auto I1_el = trace(sig_el);
+ const auto J2_el = max((s_el | s_el) / 2., local_zero_tolerance);
+ const auto J3_el = det(s_el);
+ const auto arg = min(max(-3. * sqrt3 * J3_el / (2. * J2_el * sqrt(J2_el)),
+ -1. + local_zero_tolerance),
+ 1. - local_zero_tolerance);
+ const auto lode_el = 1. / 3. * asin(arg);
+ auto K = 0.0;
+ if (fabs(lode_el) < lodeT)
+ {
+ K = cos(lode_el) - isqrt3 * sin_phi * sin(lode_el);
+ }
+ else
+ {
+ const auto sign = min(
+ max(lode_el / max(abs(lode_el), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_phi);
+ const auto B = 1 / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_phi * cos_lodeT);
+ K = A - B * arg;
+ }
+ const auto sMC =
+ I1_el / 3. * sin_phi + sqrt(J2_el * K * K + a * a * sin_phi * sin_phi);
+
+ // Maximum shear stress
+ const double t_s_el =
+ (sig_el.size() == 6)
+ ? sqrt((sig_el[3] * sig_el[3] + sig_el[4] * sig_el[4]) / 2.)
+ : sqrt((sig_el[3] * sig_el[3]) / 2);
+
+ // weak plane Coulomb
+ F[0] = t_s_el - c_wp + sig_el[0] * tan_phi_wp > 0.;
+
+ // matrix Mohr Coulomb
+ F[1] = sMC - c * cos_phi > 0.;
+ np = npwp = Stensor(real(0));
+}
+
+@Integrator
+{
+ if (F[0] || F[1])
+ {
+ const auto id = Stensor::Id();
+ const auto id4 = Stensor4::Id();
+ const auto Pdev = Stensor4::K();
+
+ if (F[0])
+ {
+ const double t_s =
+ (sig.size() == 6)
+ ? sqrt((sig[3] * sig[3] + sig[4] * sig[4]) / 2)
+ : sqrt((sig[3] * sig[3]) / 2);
+ // yield function value
+ const auto Fy = t_s - c_wp + sig[0] * tan_phi_wp;
+ if (Fy > 1e-4 * D(0, 0))
+ {
+ return false;
+ }
+
+ // flow direction and yield function gradient
+ auto n = id;
+ n *= 0.;
+ n[0] = tan_psi_wp; // x
+
+ auto nF = id;
+ nF *= 0.;
+ nF[0] = tan_phi_wp; // x
+
+ auto dn_dsig = id4;
+ dn_dsig *= 0.;
+
+ if (t_s > local_zero_tolerance)
+ {
+ n[3] = nF[3] = sig[3] / (t_s * 2.); // xy
+ const double t_s_cube_over_four = t_s * t_s * t_s * 4.;
+ dn_dsig[21] = dn_dsig[28] = 1. / (t_s * 2.); //(3,3) and (4,4)
+ dn_dsig[21] -= sig[3] * sig[3] / t_s_cube_over_four; //(3,3)
+ dn_dsig[28] -= sig[4] * sig[4] / t_s_cube_over_four; //(4,4)
+ // xz direction
+ if (sig.size() == 6)
+ n[4] = nF[4] = sig[4] / (t_s * 2); // xz
+ dn_dsig[22] = dn_dsig[24] =
+ -sig[3] * sig[4] / t_s_cube_over_four; //(3,4), (4,3)
+ }
+
+ if (this->iter > 30 &&
+ abs(n | npwp) < sqrt(n | n) * sqrt(npwp | npwp) * 0.99)
+ {
+ return false;
+ }
+
+ // residuals
+ feel += dlamWP * n;
+ flamWP = Fy / D(0, 0);
+ // Jacobian
+ dfeel_ddeel += theta * dlamWP * (dn_dsig * D);
+ dfeel_ddlamWP = n;
+ dflamWP_ddlamWP = strain(0.);
+ dflamWP_ddeel = theta * (nF | D) / D(0, 0);
+
+ npwp = n;
+ }
+
+ if (F[1])
+ {
+ constexpr auto sqrt3 = Cste<real>::sqrt3;
+ constexpr auto isqrt3 = Cste<real>::isqrt3;
+ constexpr auto id = Stensor::Id();
+ const auto s = deviator(sig);
+ const auto I1 = trace(sig);
+ const auto J2 = max((s | s) / 2., local_zero_tolerance);
+ const auto J3 =
+ real(det(s) < 0. ? min(det(s), -local_zero_tolerance)
+ : max(det(s), local_zero_tolerance));
+ const auto arg = min(max(-3. * sqrt3 * J3 / (2. * J2 * sqrt(J2)),
+ -1. + local_zero_tolerance),
+ 1. - local_zero_tolerance);
+ const auto lode = 1. / 3. * asin(arg);
+ const auto cos_lode = cos(lode);
+ const auto sin_lode = sin(lode);
+ const auto cos_3_lode = cos(3. * lode);
+ const auto sin_3_lode = arg;
+ const auto tan_3_lode = tan(3. * lode);
+ auto K = 0.;
+ auto dK_dlode = 1.;
+ if (fabs(lode) < lodeT)
+ {
+ K = cos_lode - isqrt3 * sin_phi * sin_lode;
+ dK_dlode = -sin_lode - isqrt3 * sin_phi * cos_lode;
+ }
+ else
+ {
+ const auto sign = min(
+ max(lode / max(abs(lode), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_phi);
+ const auto B =
+ 1. / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_phi * cos_lodeT);
+ K = A - B * sin_3_lode;
+ dK_dlode = -3. * B * cos_3_lode;
+ }
+ auto KG = 0.0; // move into a function to avoid code duplication
+ auto dKG_dlode = 1.;
+ auto dKG_ddlode = 1.;
+ if (fabs(lode) < lodeT)
+ {
+ KG = cos_lode - isqrt3 * sin_psi * sin_lode;
+ dKG_dlode = -sin_lode - isqrt3 * sin_psi * cos_lode;
+ dKG_ddlode = -cos_lode + isqrt3 * sin_psi * sin_lode;
+ }
+ else
+ {
+ const auto sign = min(
+ max(lode / max(fabs(lode), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_psi);
+ const auto B =
+ 1. / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_psi * cos_lodeT);
+ KG = A - B * sin_3_lode;
+ dKG_dlode = -3. * B * cos_3_lode;
+ dKG_ddlode = 9. * B * sin_3_lode;
+ }
+
+ // flow direction
+ const auto dev_s_squared = computeJ3Derivative(
+ sig); // replaces dev_s_squared = deviator(square(s));
+ const auto dG_dI1 = sin_psi / 3.;
+ const auto root =
+ max(sqrt(J2 * KG * KG +
+ a * a * tan(phi) * tan(phi) * cos(psi) * cos(psi)),
+ local_zero_tolerance);
+ const auto dG_dJ2 =
+ KG / (2. * root) * (KG - tan_3_lode * dKG_dlode);
+ const auto dG_dJ3 =
+ J2 * KG * tan_3_lode / (3. * J3 * root) * dKG_dlode;
+ const auto n =
+ eval(dG_dI1 * id + dG_dJ2 * s + dG_dJ3 * dev_s_squared);
+ if (this->iter > 30 &&
+ abs(n | np) < sqrt(n | n) * sqrt(np | np) * 0.99)
+ {
+ return false;
+ }
+ // yield function
+ const auto rootF = max(sqrt(J2 * K * K + a * a * sin_phi * sin_phi),
+ local_zero_tolerance);
+ const auto Fy = I1 * sin_phi / 3 + rootF - c * cos_phi;
+ // return if overshoot of yield surface
+ if (Fy > 1e-4 * D(0, 0))
+ {
+ return false;
+ }
+
+ // yield function derivative for Jacobian
+ const auto dF_dI1 = sin_phi / 3.;
+ const auto dF_dJ2 = K / (2. * rootF) * (K - tan_3_lode * dK_dlode);
+ const auto dF_dJ3 =
+ J2 * K * tan_3_lode / (3. * J3 * rootF) * dK_dlode;
+ const auto nF =
+ eval(dF_dI1 * id + dF_dJ2 * s + dF_dJ3 * dev_s_squared);
+
+ const auto dG_dlode = KG * J2 / (root)*dKG_dlode;
+ const auto dG_ddlode =
+ J2 / root *
+ (dKG_dlode * dKG_dlode * (1. - J2 * KG * KG / (root * root)) +
+ KG * dKG_ddlode);
+ const auto dG_ddlodeJ2 =
+ KG / root * dKG_dlode * (1. - J2 * KG * KG / (2 * root * root));
+ const auto dG_ddJ2 =
+ -KG * KG * KG * KG / (4. * root * root * root) +
+ dG_dlode * tan_3_lode / (2 * J2 * J2) -
+ tan_3_lode / (2 * J2) *
+ (2 * dG_ddlodeJ2 - tan_3_lode / (2 * J2) * dG_ddlode -
+ 3 / (2 * J2 * cos_3_lode * cos_3_lode) * dG_dlode);
+ const auto dG_ddJ3 =
+ -tan_3_lode / (3 * J3 * J3) * dG_dlode +
+ tan_3_lode / (3 * J3) *
+ (dG_ddlode * tan_3_lode / (3 * J3) +
+ dG_dlode * 1. / (J3 * cos_3_lode * cos_3_lode));
+ const auto dG_ddJ2J3 =
+ dG_ddlodeJ2 * tan_3_lode / (3 * J3) -
+ tan_3_lode / (2 * J2) *
+ (dG_ddlode * tan_3_lode / (3 * J3) +
+ dG_dlode * 1. / (J3 * cos_3_lode * cos_3_lode));
+
+ // elasticity
+ feel += dlam * n;
+ dfeel_ddeel +=
+ theta * dlam *
+ (dG_dJ2 * Pdev + dG_dJ3 * computeJ3SecondDerivative(sig) +
+ dG_ddJ2 * (s ^ s) + dG_ddJ3 * (dev_s_squared ^ dev_s_squared) +
+ dG_ddJ2J3 * ((dev_s_squared ^ s) + (s ^ dev_s_squared))) *
+ D;
+ dfeel_ddlam = n;
+ // plasticity
+ flam = Fy / D(0, 0);
+ dflam_ddlam = strain(0);
+ dflam_ddeel = theta * (nF | D) / D(0, 0);
+ np = n;
+ }
+ }
+}
+
+@AdditionalConvergenceChecks
+{
+ if (converged)
+ {
+ if (F[0])
+ {
+ if (dlamWP < 0)
+ {
+ // deactivating weak plane
+ converged = F[0] = false;
+ }
+ }
+ else
+ {
+ const double t_s =
+ (sig.size() == 6)
+ ? sqrt((sig[3] * sig[3] + sig[4] * sig[4]) / 2)
+ : sqrt((sig[3] * sig[3]) / 2);
+ // yield function value
+ const auto Fy = t_s - c_wp + sig[0] * tan_phi_wp;
+ if (Fy > 0)
+ {
+ converged = false;
+ F[0] = true;
+ }
+ }
+ if (F[1])
+ {
+ if (dlam < 0)
+ {
+ // deactivating matrix
+ converged = F[1] = false;
+ }
+ }
+ else
+ {
+ constexpr auto sqrt3 = Cste<real>::sqrt3;
+ constexpr auto isqrt3 = Cste<real>::isqrt3;
+ const auto s = deviator(sig);
+ const auto I1 = trace(sig);
+ const auto J2 = max((s | s) / 2., local_zero_tolerance);
+ const auto J3 =
+ real(det(s) < 0. ? min(det(s), -local_zero_tolerance)
+ : max(det(s), local_zero_tolerance));
+ const auto arg = min(max(-3. * sqrt3 * J3 / (2. * J2 * sqrt(J2)),
+ -1. + local_zero_tolerance),
+ 1. - local_zero_tolerance);
+ const auto lode = 1. / 3. * asin(arg);
+ const auto cos_lode = cos(lode);
+ const auto sin_lode = sin(lode);
+ const auto sin_3_lode = arg;
+ auto K = 0.;
+ if (fabs(lode) < lodeT)
+ {
+ K = cos_lode - isqrt3 * sin_phi * sin_lode;
+ }
+ else
+ {
+ const auto sign = min(
+ max(lode / max(abs(lode), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_phi);
+ const auto B =
+ 1. / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_phi * cos_lodeT);
+ K = A - B * sin_3_lode;
+ }
+ // yield function
+ const auto rootF = max(sqrt(J2 * K * K + a * a * sin_phi * sin_phi),
+ local_zero_tolerance);
+ const auto Fy = I1 * sin_phi / 3 + rootF - c * cos_phi;
+ if (Fy > 0)
+ {
+ converged = false;
+ F[1] = true;
+ }
+ }
+ }
+}
diff --git a/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanUBI.mtest b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanUBI.mtest
new file mode 100644
index 0000000000000000000000000000000000000000..3a9ca1e0816e100bf5c895719788442a8b1d92a6
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanUBI.mtest
@@ -0,0 +1,40 @@
+@Author Thomas Nagel;
+@Date 19/02/2019;
+@Description{
+
+};
+
+
+@Real 'theta' -2.697327025809418;
+@Real 'em' 1e-2;
+@Real 'tmax' 1;
+@Real 'c' 'cos(theta)';
+@Real 's' 'sin(theta)';
+
+@MaximumNumberOfSubSteps 5; // fails without substepping (residual ping-pong)
+@Behaviour<generic> 'src/libBehaviour.so' 'MohrCoulombAbboSloanUBI';
+
+@ImposedStrain 'EXX' {0 : 0, 'tmax' : 'em* c'};
+@ImposedStrain 'EYY' {0 : 0, 'tmax' : 'em* s'};
+@ImposedStrain 'EXY' {0 : 0, 0.5 : 0.0, 1.0 : 'em'};
+@ImposedStrain 'EXZ' {0 : 0, 0.5 : 0.0, 1.0 : '2*em'};
+
+//@NonLinearConstraint<Stress> 'SXX+SYY+SZZ';
+//@NonLinearConstraint<Stress> 'SXY';
+//@NonLinearConstraint<Stress> 'SXZ';
+//@NonLinearConstraint<Stress> 'SYZ';
+
+@ExternalStateVariable "Temperature" 293.15;
+
+@Times{0, 1 in 100};
+
+@MaterialProperty<constant> 'YoungModulus' 150.e3;
+@MaterialProperty<constant> 'PoissonRatio' 0.3;
+@MaterialProperty<constant> 'PlaneCohesion' 2.e1;
+@MaterialProperty<constant> 'PlaneFrictionAngle' 10.;
+@MaterialProperty<constant> 'PlaneDilatancyAngle' 10.;
+@MaterialProperty<constant> 'Cohesion' 3.e1;
+@MaterialProperty<constant> 'FrictionAngle' 30.;
+@MaterialProperty<constant> 'DilatancyAngle' 10.;
+@MaterialProperty<constant> 'TransitionAngle' 29.;
+@MaterialProperty<constant> 'TensionCutOffParameter' 1.e1;
diff --git a/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanUBIOrtho.mfront b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanUBIOrtho.mfront
new file mode 100644
index 0000000000000000000000000000000000000000..1d418bacd8003cdaa75d5a4b1be68eaaa0dab86e
--- /dev/null
+++ b/MaterialLib/SolidModels/MFront/MohrCoulombAbboSloanUBIOrtho.mfront
@@ -0,0 +1,442 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+@DSL Implicit;
+@Behaviour MohrCoulombAbboSloanUBIOrtho;
+@Author Thomas Nagel;
+@Date 21 / 03 / 2020;
+@Description{
+Ubiquitous joint model where the matrix fails according to a Mohr-Coulomb
+criterion(with tension-cutoff) and the embedded weakness planes follow
+a Coulomb behaviour(so far without tension cut-off). Both models are
+non-associated.
+The fracture normal is the x-direction(1, 0, 0).
+Any rotations to be done by the FE programme.
+The elastic stiffness tensor is orthotropic.
+};
+
+@Algorithm NewtonRaphson;
+@MaximumNumberOfIterations 100;
+
+@OrthotropicBehaviour<Pipe>;
+@Brick StandardElasticity;
+
+@Epsilon 1.e-14;
+@Theta 1.0;
+@Parameter local_zero_tolerance = 1.e-14;
+@Parameter pi = 3.14159265359;
+
+@ModellingHypotheses{".+"};
+@RequireStiffnessTensor<UnAltered>;
+
+// weakness plane cohesion
+@MaterialProperty real c_wp;
+c_wp.setEntryName("PlaneCohesion");
+// weakness plane friction angle
+@MaterialProperty real phi_wp;
+phi_wp.setEntryName("PlaneFrictionAngle");
+// weakness plane dilatancy angle
+@MaterialProperty real psi_wp;
+psi_wp.setEntryName("PlaneDilatancyAngle");
+
+@MaterialProperty stress c;
+c.setEntryName("Cohesion");
+@MaterialProperty real phi;
+phi.setEntryName("FrictionAngle");
+@MaterialProperty real psi;
+psi.setEntryName("DilatancyAngle");
+@MaterialProperty real lodeT;
+lodeT.setEntryName("TransitionAngle");
+@MaterialProperty stress a;
+a.setEntryName("TensionCutOffParameter");
+
+@LocalVariable Stensor np;
+@LocalVariable Stensor npwp;
+
+// Not array because of OGS output
+@StateVariable strain lamWP;
+lamWP.setEntryName("EquivalentPlasticStrainWP");
+@StateVariable strain lam;
+lam.setEntryName("EquivalentPlasticStrainMatrix");
+
+@LocalVariable bool F[2];
+@LocalVariable real tan_phi_wp;
+@LocalVariable real tan_psi_wp;
+@LocalVariable real sin_psi;
+@LocalVariable real sin_phi;
+@LocalVariable real cos_phi;
+@LocalVariable real cos_lodeT;
+@LocalVariable real sin_lodeT;
+@LocalVariable real tan_lodeT;
+@LocalVariable real cos_3_lodeT;
+@LocalVariable real sin_3_lodeT;
+@LocalVariable real tan_3_lodeT;
+
+@InitLocalVariables
+{
+ // tan_phi_wp after conversion to rad
+ tan_phi_wp = tan(phi_wp * pi / 180.);
+ // tan_psi_wp after conversion to rad
+ tan_psi_wp = tan(psi_wp * pi / 180.);
+
+ constexpr auto sqrt3 = Cste<real>::sqrt3;
+ constexpr auto isqrt3 = Cste<real>::isqrt3;
+ // conversion to rad
+ phi *= pi / 180.;
+ psi *= pi / 180.;
+ lodeT *= pi / 180.;
+ sin_psi = sin(psi);
+ cos_phi = cos(phi);
+ sin_phi = sin(phi);
+ sin_lodeT = sin(lodeT);
+ cos_lodeT = cos(lodeT);
+ tan_lodeT = tan(lodeT);
+ cos_3_lodeT = cos(3. * lodeT);
+ sin_3_lodeT = sin(3. * lodeT);
+ tan_3_lodeT = tan(3. * lodeT);
+
+ // Compute initial elastic strain
+ const auto S = invert(D);
+ eel = S * sig;
+
+ const auto sig_el = computeElasticPrediction();
+
+ const auto s_el = deviator(sig_el);
+ const auto I1_el = trace(sig_el);
+ const auto J2_el = max((s_el | s_el) / 2., local_zero_tolerance);
+ const auto J3_el = det(s_el);
+ const auto arg = min(max(-3. * sqrt3 * J3_el / (2. * J2_el * sqrt(J2_el)),
+ -1. + local_zero_tolerance),
+ 1. - local_zero_tolerance);
+ const auto lode_el = 1. / 3. * asin(arg);
+ auto K = 0.0;
+ if (fabs(lode_el) < lodeT)
+ {
+ K = cos(lode_el) - isqrt3 * sin_phi * sin(lode_el);
+ }
+ else
+ {
+ const auto sign = min(
+ max(lode_el / max(abs(lode_el), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_phi);
+ const auto B = 1 / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_phi * cos_lodeT);
+ K = A - B * arg;
+ }
+ const auto sMC =
+ I1_el / 3. * sin_phi + sqrt(J2_el * K * K + a * a * sin_phi * sin_phi);
+
+ // Maximum shear stress
+ const double t_s_el =
+ (sig_el.size() == 6)
+ ? sqrt((sig_el[3] * sig_el[3] + sig_el[4] * sig_el[4]) / 2.)
+ : sqrt((sig_el[3] * sig_el[3]) / 2);
+
+ // weak plane Coulomb
+ F[0] = t_s_el - c_wp + sig_el[0] * tan_phi_wp > 0.;
+
+ // matrix Mohr Coulomb
+ F[1] = sMC - c * cos_phi > 0.;
+ np = npwp = Stensor(real(0));
+}
+
+@Integrator
+{
+ if (F[0] || F[1])
+ {
+ const auto id = Stensor::Id();
+ const auto id4 = Stensor4::Id();
+ const auto Pdev = Stensor4::K();
+
+ if (F[0])
+ {
+ const double t_s =
+ (sig.size() == 6)
+ ? sqrt((sig[3] * sig[3] + sig[4] * sig[4]) / 2)
+ : sqrt((sig[3] * sig[3]) / 2);
+ // yield function value
+ const auto Fy = t_s - c_wp + sig[0] * tan_phi_wp;
+ if (Fy > 1e-4 * D(0, 0))
+ {
+ return false;
+ }
+
+ // flow direction and yield function gradient
+ auto n = id;
+ n *= 0.;
+ n[0] = tan_psi_wp; // x
+
+ auto nF = id;
+ nF *= 0.;
+ nF[0] = tan_phi_wp; // x
+
+ auto dn_dsig = id4;
+ dn_dsig *= 0.;
+
+ if (t_s > local_zero_tolerance)
+ {
+ n[3] = nF[3] = sig[3] / (t_s * 2.); // xy
+ const double t_s_cube_over_four = t_s * t_s * t_s * 4.;
+ dn_dsig[21] = dn_dsig[28] = 1. / (t_s * 2.); //(3,3) and (4,4)
+ dn_dsig[21] -= sig[3] * sig[3] / t_s_cube_over_four; //(3,3)
+ dn_dsig[28] -= sig[4] * sig[4] / t_s_cube_over_four; //(4,4)
+ // xz direction
+ if (sig.size() == 6)
+ n[4] = nF[4] = sig[4] / (t_s * 2); // xz
+ dn_dsig[22] = dn_dsig[24] =
+ -sig[3] * sig[4] / t_s_cube_over_four; //(3,4), (4,3)
+ }
+
+ if (this->iter > 30 &&
+ abs(n | npwp) < sqrt(n | n) * sqrt(npwp | npwp) * 0.99)
+ {
+ return false;
+ }
+
+ // residuals
+ feel += dlamWP * n;
+ flamWP = Fy / D(0, 0);
+ // Jacobian
+ dfeel_ddeel += theta * dlamWP * (dn_dsig * D);
+ dfeel_ddlamWP = n;
+ dflamWP_ddlamWP = strain(0.);
+ dflamWP_ddeel = theta * (nF | D) / D(0, 0);
+
+ npwp = n;
+ }
+
+ if (F[1])
+ {
+ constexpr auto sqrt3 = Cste<real>::sqrt3;
+ constexpr auto isqrt3 = Cste<real>::isqrt3;
+ constexpr auto id = Stensor::Id();
+ const auto s = deviator(sig);
+ const auto I1 = trace(sig);
+ const auto J2 = max((s | s) / 2., local_zero_tolerance);
+ const auto J3 =
+ real(det(s) < 0. ? min(det(s), -local_zero_tolerance)
+ : max(det(s), local_zero_tolerance));
+ const auto arg = min(max(-3. * sqrt3 * J3 / (2. * J2 * sqrt(J2)),
+ -1. + local_zero_tolerance),
+ 1. - local_zero_tolerance);
+ const auto lode = 1. / 3. * asin(arg);
+ const auto cos_lode = cos(lode);
+ const auto sin_lode = sin(lode);
+ const auto cos_3_lode = cos(3. * lode);
+ const auto sin_3_lode = arg;
+ const auto tan_3_lode = tan(3. * lode);
+ auto K = 0.;
+ auto dK_dlode = 1.;
+ if (fabs(lode) < lodeT)
+ {
+ K = cos_lode - isqrt3 * sin_phi * sin_lode;
+ dK_dlode = -sin_lode - isqrt3 * sin_phi * cos_lode;
+ }
+ else
+ {
+ const auto sign = min(
+ max(lode / max(abs(lode), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_phi);
+ const auto B =
+ 1. / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_phi * cos_lodeT);
+ K = A - B * sin_3_lode;
+ dK_dlode = -3. * B * cos_3_lode;
+ }
+ auto KG = 0.0; // move into a function to avoid code duplication
+ auto dKG_dlode = 1.;
+ auto dKG_ddlode = 1.;
+ if (fabs(lode) < lodeT)
+ {
+ KG = cos_lode - isqrt3 * sin_psi * sin_lode;
+ dKG_dlode = -sin_lode - isqrt3 * sin_psi * cos_lode;
+ dKG_ddlode = -cos_lode + isqrt3 * sin_psi * sin_lode;
+ }
+ else
+ {
+ const auto sign = min(
+ max(lode / max(fabs(lode), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_psi);
+ const auto B =
+ 1. / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_psi * cos_lodeT);
+ KG = A - B * sin_3_lode;
+ dKG_dlode = -3. * B * cos_3_lode;
+ dKG_ddlode = 9. * B * sin_3_lode;
+ }
+
+ // flow direction
+ const auto dev_s_squared = computeJ3Derivative(
+ sig); // replaces dev_s_squared = deviator(square(s));
+ const auto dG_dI1 = sin_psi / 3.;
+ const auto root =
+ max(sqrt(J2 * KG * KG +
+ a * a * tan(phi) * tan(phi) * cos(psi) * cos(psi)),
+ local_zero_tolerance);
+ const auto dG_dJ2 =
+ KG / (2. * root) * (KG - tan_3_lode * dKG_dlode);
+ const auto dG_dJ3 =
+ J2 * KG * tan_3_lode / (3. * J3 * root) * dKG_dlode;
+ const auto n =
+ eval(dG_dI1 * id + dG_dJ2 * s + dG_dJ3 * dev_s_squared);
+ if (this->iter > 30 &&
+ abs(n | np) < sqrt(n | n) * sqrt(np | np) * 0.99)
+ {
+ return false;
+ }
+ // yield function
+ const auto rootF = max(sqrt(J2 * K * K + a * a * sin_phi * sin_phi),
+ local_zero_tolerance);
+ const auto Fy = I1 * sin_phi / 3 + rootF - c * cos_phi;
+ // return if overshoot of yield surface
+ if (Fy > 1e-4 * D(0, 0))
+ {
+ return false;
+ }
+
+ // yield function derivative for Jacobian
+ const auto dF_dI1 = sin_phi / 3.;
+ const auto dF_dJ2 = K / (2. * rootF) * (K - tan_3_lode * dK_dlode);
+ const auto dF_dJ3 =
+ J2 * K * tan_3_lode / (3. * J3 * rootF) * dK_dlode;
+ const auto nF =
+ eval(dF_dI1 * id + dF_dJ2 * s + dF_dJ3 * dev_s_squared);
+
+ const auto dG_dlode = KG * J2 / (root)*dKG_dlode;
+ const auto dG_ddlode =
+ J2 / root *
+ (dKG_dlode * dKG_dlode * (1. - J2 * KG * KG / (root * root)) +
+ KG * dKG_ddlode);
+ const auto dG_ddlodeJ2 =
+ KG / root * dKG_dlode * (1. - J2 * KG * KG / (2 * root * root));
+ const auto dG_ddJ2 =
+ -KG * KG * KG * KG / (4. * root * root * root) +
+ dG_dlode * tan_3_lode / (2 * J2 * J2) -
+ tan_3_lode / (2 * J2) *
+ (2 * dG_ddlodeJ2 - tan_3_lode / (2 * J2) * dG_ddlode -
+ 3 / (2 * J2 * cos_3_lode * cos_3_lode) * dG_dlode);
+ const auto dG_ddJ3 =
+ -tan_3_lode / (3 * J3 * J3) * dG_dlode +
+ tan_3_lode / (3 * J3) *
+ (dG_ddlode * tan_3_lode / (3 * J3) +
+ dG_dlode * 1. / (J3 * cos_3_lode * cos_3_lode));
+ const auto dG_ddJ2J3 =
+ dG_ddlodeJ2 * tan_3_lode / (3 * J3) -
+ tan_3_lode / (2 * J2) *
+ (dG_ddlode * tan_3_lode / (3 * J3) +
+ dG_dlode * 1. / (J3 * cos_3_lode * cos_3_lode));
+
+ // elasticity
+ feel += dlam * n;
+ dfeel_ddeel +=
+ theta * dlam *
+ (dG_dJ2 * Pdev + dG_dJ3 * computeJ3SecondDerivative(sig) +
+ dG_ddJ2 * (s ^ s) + dG_ddJ3 * (dev_s_squared ^ dev_s_squared) +
+ dG_ddJ2J3 * ((dev_s_squared ^ s) + (s ^ dev_s_squared))) *
+ D;
+ dfeel_ddlam = n;
+ // plasticity
+ flam = Fy / D(0, 0);
+ dflam_ddlam = strain(0);
+ dflam_ddeel = theta * (nF | D) / D(0, 0);
+ np = n;
+ }
+ }
+}
+
+@AdditionalConvergenceChecks
+{
+ if (converged)
+ {
+ if (F[0])
+ {
+ if (dlamWP < 0)
+ {
+ // deactivating weak plane
+ converged = F[0] = false;
+ }
+ }
+ else
+ {
+ const double t_s =
+ (sig.size() == 6)
+ ? sqrt((sig[3] * sig[3] + sig[4] * sig[4]) / 2)
+ : sqrt((sig[3] * sig[3]) / 2);
+ // yield function value
+ const auto Fy = t_s - c_wp + sig[0] * tan_phi_wp;
+ if (Fy > 0)
+ {
+ converged = false;
+ F[0] = true;
+ }
+ }
+ if (F[1])
+ {
+ if (dlam < 0)
+ {
+ // deactivating matrix
+ converged = F[1] = false;
+ }
+ }
+ else
+ {
+ constexpr auto sqrt3 = Cste<real>::sqrt3;
+ constexpr auto isqrt3 = Cste<real>::isqrt3;
+ const auto s = deviator(sig);
+ const auto I1 = trace(sig);
+ const auto J2 = max((s | s) / 2., local_zero_tolerance);
+ const auto J3 =
+ real(det(s) < 0. ? min(det(s), -local_zero_tolerance)
+ : max(det(s), local_zero_tolerance));
+ const auto arg = min(max(-3. * sqrt3 * J3 / (2. * J2 * sqrt(J2)),
+ -1. + local_zero_tolerance),
+ 1. - local_zero_tolerance);
+ const auto lode = 1. / 3. * asin(arg);
+ const auto cos_lode = cos(lode);
+ const auto sin_lode = sin(lode);
+ const auto sin_3_lode = arg;
+ auto K = 0.;
+ if (fabs(lode) < lodeT)
+ {
+ K = cos_lode - isqrt3 * sin_phi * sin_lode;
+ }
+ else
+ {
+ const auto sign = min(
+ max(lode / max(abs(lode), local_zero_tolerance), -1.), 1.);
+ const auto A =
+ 1. / 3. * cos_lodeT *
+ (3. + tan_lodeT * tan_3_lodeT +
+ isqrt3 * sign * (tan_3_lodeT - 3. * tan_lodeT) * sin_phi);
+ const auto B =
+ 1. / (3. * cos_3_lodeT) *
+ (sign * sin_lodeT + isqrt3 * sin_phi * cos_lodeT);
+ K = A - B * sin_3_lode;
+ }
+ // yield function
+ const auto rootF = max(sqrt(J2 * K * K + a * a * sin_phi * sin_phi),
+ local_zero_tolerance);
+ const auto Fy = I1 * sin_phi / 3 + rootF - c * cos_phi;
+ if (Fy > 0)
+ {
+ converged = false;
+ F[1] = true;
+ }
+ }
+ }
+}
diff --git a/ProcessLib/SmallDeformation/Tests.cmake b/ProcessLib/SmallDeformation/Tests.cmake
index a804d38feb6f380f150d5d64368a7507e83e187a..195ce4f9d78ceb13ef9e3ffd4347d4b41fc4a894 100644
--- a/ProcessLib/SmallDeformation/Tests.cmake
+++ b/ProcessLib/SmallDeformation/Tests.cmake
@@ -52,6 +52,9 @@ endif()
if (OGS_USE_MFRONT)
OgsTest(PROJECTFILE Mechanics/MohrCoulombAbboSloan/load_test_mc.prj)
+ OgsTest(PROJECTFILE Mechanics/MohrCoulombAbboSloanAnisotropic/triax_1e0_47.prj)
+ OgsTest(PROJECTFILE Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ortho_1e0_47.prj)
+ OgsTest(PROJECTFILE Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_1e0_47.prj)
#TODO (naumov) enable when output file format can be specified
#OgsTest(PROJECTFILE Mechanics/MohrCoulombAbboSloan/oedometer.prj RUNTIME 80)
OgsTest(PROJECTFILE Mechanics/Linear/MFront/cube_1e0_orthotropic_xyz.prj)
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/run_script.py b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/run_script.py
new file mode 100644
index 0000000000000000000000000000000000000000..8fc91d87d11b3889483e70eb9fce285b1630964d
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/run_script.py
@@ -0,0 +1,143 @@
+#!/ usr / bin / env python
+#coding : utf - 8
+
+#
+#
+#Run skript for Mohr Coulomb with anisotropy via stress scaling or \
+ embedded weakness planes.
+#Comments to:
+#
+#* Prof.Dr.Thomas Nagel
+#Chair of Soil Mechanics and Foundation Engineering
+#Geotechnical Institute
+#Technische Universität Bergakademie Freiberg.*
+#
+#https: // tu-freiberg.de/en/fakultaet3/gt/soilmechanics
+#
+#Tests and parameter sets based on:
+#
+#Malcom, M.A.M.(2018).Analysis of underground excavations in argillaceous hard \
+ soils - \
+ weak rocks.Technical University of Catalonia.
+#
+#Ismael, M., &Konietzky, \
+ H.(2017) \
+ .Integration of Elastic Stiffness Anisotropy into Ubiquitous Joint \
+ Model.Procedia Engineering, \
+ 191, 1032–1039. https: // doi.org/10.1016/j.proeng.2017.05.276
+
+#In[1]:
+
+import numpy as np import os as os import pyvista as pv import matplotlib.pyplot as plt
+
+#In[2]:
+
+#Some plot settings
+ plt.style.use('seaborn-deep') plt.rcParams['lines.linewidth'] = 2.0 plt.rcParams['lines.color'] = 'black' plt.rcParams['legend.frameon'] = True plt.rcParams['font.family'] = 'serif' plt.rcParams['legend.fontsize'] = 14 plt.rcParams['font.size'] = 14 plt.rcParams['axes.spines.right'] = False plt.rcParams['axes.spines.top'] = False plt.rcParams['axes.spines.left'] = True plt.rcParams['axes.spines.bottom'] = True plt.rcParams['axes.axisbelow'] = True plt.rcParams['figure.figsize'] =(12, 6)
+
+#In[3]:
+
+ cmd = 'rm triax*.pvd triax_*pcs_0*.vtu triax_*1e0_*.prj' os.system(cmd)
+
+#In[4]:
+
+ beta = np.linspace(0, np.pi / 2, 28)
+
+#In[5]:
+
+ for i in beta : print("Generating input for %i °" %(int(np.round(i * 180 / np.pi, 0)))) filename = "triax_1e0_" + str(int(np.round(i * 180. / np.pi, 0))) + ".prj" cmd = "sed 's/<values>1 0/<values>" + str(np.cos(i)) + " " + str(np.sin(i)) + "/g' triax_original.prj > " + filename os.system(cmd) cmd = "sed 's/<values>0 1/<values>" + str(- np.sin(i)) + " " + str(np.cos(i)) + "/g' " + filename + " -i" os.system(cmd) cmd = "sed 's/<prefix>triax_0/<prefix>triax_" + str(int(np.round(i * 180. / np.pi, 0))) + "/g' " + filename + " -i" os.system(cmd) print("Running simulation for %i °" %(int(np.round(i * 180 / np.pi, 0)))) cmd = "~/ogs_release/bin/ogs " + filename os.system(cmd)
+
+#In[6]:
+
+ UCS = beta * 0. for n, i in enumerate(beta) :
+#star, as ts value unknown due to adaptive time stepping
+ a = 'triax_' + str(int(np.round(i * 180. / np.pi, 0))) + '_pcs_0_ts_*_t_2.000000.vtu' filename = get_ipython().getoutput('ls $a') data = pv.read(filename[0]) sigma = data.get_array('sigma') UCS[n] = - sigma[0][1]
+
+#In[7]:
+
+ fig, ax = plt.subplots() ax.plot(90 - beta * 180 / np.pi, UCS / 1e3, marker = 'd') ax.set_xlabel('Angle of fracture plane / °') ax.set_ylabel('UCS / kPa') ax.set_ylim(2, 10) ax.set_xlim(0, 90);
+ax.set_title('isotropic elasticity')
+fig.savefig('orientation_dependent_strength.pdf')
+
+#In[8]:
+
+
+for i in beta:
+ print("Generating input for %i °" %(int(np.round(i*180/np.pi,0))))
+ filename = "triax_ortho_1e0_"+str(int(np.round(i*180./np.pi,0)))+".prj"
+ cmd = "sed 's/<values>1 0/<values>" + str(np.cos(i)) + " " + str(np.sin(i))+ "/g' triax_ORTHO_original.prj > " + filename
+ os.system(cmd)
+ cmd = "sed 's/<values>0 1/<values>" + str(-np.sin(i)) + " " + str(np.cos(i))+ "/g' " + filename + " -i"
+ os.system(cmd)
+ cmd = "sed 's/<prefix>triax_ortho_0/<prefix>triax_ortho_" +str(int(np.round(i*180./np.pi,0))) + "/g' " + filename + " -i"
+ os.system(cmd)
+ print("Running simulation for %i °" %(int(np.round(i*180/np.pi,0))))
+ cmd = "~/ogs_release/bin/ogs " + filename
+ os.system(cmd)
+
+#In[9]:
+
+
+UCS = beta * 0.
+for n,i in enumerate(beta):
+#star, as ts value unknown due to adaptive time stepping
+ a = 'triax_ortho_'+str(int(np.round(i*180./np.pi,0)))+'_pcs_0_ts_*_t_2.000000.vtu'
+ filename = get_ipython().getoutput('ls $a')
+ data = pv.read(filename[0])
+ sigma = data.get_array('sigma')
+ UCS[n] = -sigma[0][1]
+
+#In[10]:
+
+
+fig, ax = plt.subplots()
+ax.plot(90 - beta*180/np.pi,UCS/1e6,marker='d')
+ax.set_xlabel('Angle of fracture plane / °')
+ax.set_ylabel('UCS / MPa')
+ax.set_ylim(0,120)
+ax.set_xlim(0,90);
+ax.set_title('transversely isotropic elasticity')
+fig.savefig('orientation_dependent_strength_ortho.pdf')
+
+#In[11]:
+
+
+for i in beta:
+ print("Generating input for %i °" %(int(np.round(i*180/np.pi,0))))
+ filename = "triax_aniso_1e0_"+str(int(np.round(i*180./np.pi,0)))+".prj"
+ cmd = "sed 's/<values>1 0/<values>" + str(np.cos(i)) + " " + str(np.sin(i))+ "/g' triax_aniso_original.prj > " + filename
+ os.system(cmd)
+ cmd = "sed 's/<values>0 1/<values>" + str(-np.sin(i)) + " " + str(np.cos(i))+ "/g' " + filename + " -i"
+ os.system(cmd)
+ cmd = "sed 's/<prefix>triax_aniso_0/<prefix>triax_aniso_" +str(int(np.round(i*180./np.pi,0))) + "/g' " + filename + " -i"
+ os.system(cmd)
+ print("Running simulation for %i °" %(int(np.round(i*180/np.pi,0))))
+ cmd = "~/ogs_release/bin/ogs " + filename
+ os.system(cmd)
+
+#In[12]:
+
+
+UCS = beta * 0.
+for n,i in enumerate(beta):
+#star, as ts value unknown due to adaptive time stepping
+ a = 'triax_aniso_'+str(int(np.round(i*180./np.pi,0)))+'_pcs_0_ts_*_t_2.000000.vtu'
+ filename = get_ipython().getoutput('ls $a')
+ data = pv.read(filename[0])
+ sigma = data.get_array('sigma')
+ UCS[n] = -sigma[0][1]
+
+#In[13]:
+
+
+fig, ax = plt.subplots()
+ax.plot(90 - beta*180/np.pi,UCS/1e6,marker='d')
+ax.set_xlabel('Angle of fracture plane / °')
+ax.set_ylabel('UCS / MPa')
+ax.set_ylim(0,120)
+ax.set_xlim(0,90);
+ax.set_title('transversely isotropic elasticity') fig.savefig(
+ 'orientation_dependent_strength_aniso.pdf')
+
+#In[]:
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/square_1x1.gml b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/square_1x1.gml
new file mode 100644
index 0000000000000000000000000000000000000000..4e221adcca321d6abc2bffb6cfa6da6f34366c79
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/square_1x1.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6ce1406824f519ea96f18a28e0db575a8d71dd05332b6108af721bf985e9c250
+size 944
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/square_1x1_quad_1e0.vtu b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/square_1x1_quad_1e0.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..2e6fc2e4412db4a41dd4be8bf0f0f0b3f2d94ea2
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/square_1x1_quad_1e0.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:cd87128b2c36d59cff91fd9d072f97be311ecc5119a6e607fd58e318e25b441c
+size 1578
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_1e0_47.prj b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_1e0_47.prj
new file mode 100644
index 0000000000000000000000000000000000000000..5b23ebe9ad8fbcfa7cdd214d389ea2bcd3e8f60d
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_1e0_47.prj
@@ -0,0 +1,287 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>square_1x1_quad_1e0.vtu</mesh>
+ <geometry>square_1x1.gml</geometry>
+ <processes>
+ <process>
+ <name>SD</name>
+ <type>SMALL_DEFORMATION</type>
+ <integration_order>2</integration_order>
+ <constitutive_relation>
+ <type>MFront</type>
+ <behaviour>MohrCoulombAbboSloanUBI</behaviour>
+ <material_properties>
+ <material_property name="YoungModulus" parameter="YoungModulus"/>
+ <material_property name="PoissonRatio" parameter="PoissonRatio"/>
+ <material_property name="PlaneCohesion" parameter="PlaneCohesion"/>
+ <material_property name="PlaneFrictionAngle" parameter="PlaneFrictionAngle"/>
+ <material_property name="PlaneDilatancyAngle" parameter="PlaneDilatancyAngle"/>
+ <material_property name="Cohesion" parameter="Cohesion"/>
+ <material_property name="FrictionAngle" parameter="FrictionAngle"/>
+ <material_property name="DilatancyAngle" parameter="DilatancyAngle"/>
+ <material_property name="TransitionAngle" parameter="TransitionAngle"/>
+ <material_property name="TensionCutOffParameter" parameter="TensionCutOffParameter"/>
+ </material_properties>
+ </constitutive_relation>
+ <solid_density>rho_sr</solid_density>
+ <specific_body_force>0 0</specific_body_force>
+ <process_variables>
+ <process_variable>displacement</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="ElasticStrain" output_name="ElasticStrain"/>
+ <!--secondary_variable type="static" internal_name="EquivalentPlasticStrain" output_name="EquivalentPlasticStrain"/-->
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="SD">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1e-14</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>IterationNumberBasedTimeStepping</type>
+ <t_initial>0.0</t_initial>
+ <t_end>2.0</t_end>
+ <initial_dt>0.01</initial_dt>
+ <minimum_dt>0.01</minimum_dt>
+ <maximum_dt>0.01</maximum_dt>
+ <number_iterations>1 4 10 20 </number_iterations>
+ <multiplier>1.2 1.0 0.9 0.8</multiplier>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>triax_47_pcs_{:process_id}</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <each_steps>10000</each_steps>
+ </pair>
+ </timesteps>
+ <variables>
+ <variable>displacement</variable>
+ <variable>sigma</variable>
+ </variables>
+ <suffix>_ts_{:timestep}_t_{:time}</suffix>
+ </output>
+ </time_loop>
+ <local_coordinate_system>
+ <basis_vector_0>e0</basis_vector_0>
+ <basis_vector_1>e1</basis_vector_1>
+ </local_coordinate_system>
+ <parameters>
+ <parameter>
+ <name>e0</name>
+ <type>Constant</type>
+ <values>0.6862416378687336 0.7273736415730486</values>
+ </parameter>
+ <parameter>
+ <name>e1</name>
+ <type>Constant</type>
+ <values>-0.7273736415730486 0.6862416378687336</values>
+ </parameter>
+ <parameter>
+ <name>YoungModulus</name>
+ <type>Constant</type>
+ <value>170e6</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio</name>
+ <type>Constant</type>
+ <value>.22</value>
+ </parameter>
+ <parameter>
+ <name>PlaneCohesion</name>
+ <type>Constant</type>
+ <value>1e3</value>
+ </parameter>
+ <parameter>
+ <name>PlaneFrictionAngle</name>
+ <type>Constant</type>
+ <value>30</value>
+ </parameter>
+ <parameter>
+ <name>PlaneDilatancyAngle</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>Cohesion</name>
+ <type>Constant</type>
+ <value>2.e3</value>
+ </parameter>
+ <parameter>
+ <name>FrictionAngle</name>
+ <type>Constant</type>
+ <value>40</value>
+ </parameter>
+ <parameter>
+ <name>DilatancyAngle</name>
+ <type>Constant</type>
+ <value>10</value>
+ </parameter>
+ <parameter>
+ <name>TransitionAngle</name>
+ <type>Constant</type>
+ <value>29</value>
+ </parameter>
+ <parameter>
+ <name>TensionCutOffParameter</name>
+ <type>Constant</type>
+ <value>1.0</value>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0</values>
+ </parameter>
+ <parameter>
+ <name>zero</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>axial_displacement_top</name>
+ <type>CurveScaled</type>
+ <curve>dis_loading_curve</curve>
+ <parameter>loading_value_top</parameter>
+ </parameter>
+ <parameter>
+ <name>loading_value_top</name>
+ <type>Constant</type>
+ <value>-1.e-4</value>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_left</name>
+ <type>Constant</type>
+ <values>1e0</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_left</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_left</parameter>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_right</name>
+ <type>Constant</type>
+ <values>-1e0</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_right</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_right</parameter>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>dis_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 0.0 1.0 </values>
+ </curve>
+ <curve>
+ <name>pre_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 1.0 1.0 </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>square_1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>origin</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>axial_displacement_top</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_right</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_left</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>40</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>SparseLU</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ </eigen>
+ <petsc>
+ <prefix>sd</prefix>
+ <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+ <test_definition>
+ <vtkdiff>
+ <file>triax_47_pcs_0_ts_200_t_2.000000.vtu</file>
+ <field>displacement</field>
+ <absolute_tolerance>1e-14</absolute_tolerance>
+ <relative_tolerance>1e-14</relative_tolerance>
+ </vtkdiff>
+ <vtkdiff>
+ <file>triax_47_pcs_0_ts_200_t_2.000000.vtu</file>
+ <field>sigma</field>
+ <!-- E is 1e8 -> E * u_tol = 1e-6 -->
+ <absolute_tolerance>1e-6</absolute_tolerance>
+ <relative_tolerance>1e-6</relative_tolerance>
+ </vtkdiff>
+ </test_definition>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_47_pcs_0_ts_200_t_2.000000.vtu b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_47_pcs_0_ts_200_t_2.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..985eef23e058e6e1ce863914ca7a37822f2b86e4
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_47_pcs_0_ts_200_t_2.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4f74250e2bb3a1377ce81dc1860263ac3d4fde8e75ba2c618636e5c2f57211ca
+size 5312
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ORTHO_original.prj b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ORTHO_original.prj
new file mode 100644
index 0000000000000000000000000000000000000000..c2d15d787fcdeec47e857d1366495def4a562f35
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ORTHO_original.prj
@@ -0,0 +1,314 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>square_1x1_quad_1e0.vtu</mesh>
+ <geometry>square_1x1.gml</geometry>
+ <processes>
+ <process>
+ <name>SD</name>
+ <type>SMALL_DEFORMATION</type>
+ <integration_order>2</integration_order>
+ <constitutive_relation>
+ <type>MFront</type>
+ <behaviour>MohrCoulombAbboSloanUBIOrtho</behaviour>
+ <material_properties>
+ <material_property name="YoungModulus1" parameter="YoungModulus1"/>
+ <material_property name="YoungModulus2" parameter="YoungModulus2"/>
+ <material_property name="YoungModulus3" parameter="YoungModulus3"/>
+ <material_property name="PoissonRatio12" parameter="PoissonRatio12"/>
+ <material_property name="PoissonRatio23" parameter="PoissonRatio23"/>
+ <material_property name="PoissonRatio13" parameter="PoissonRatio13"/>
+ <material_property name="ShearModulus12" parameter="ShearModulus12"/>
+ <!--material_property name="ShearModulus23" parameter="ShearModulus23"/>
+ <material_property name="ShearModulus13" parameter="ShearModulus13"/-->
+ <material_property name="PlaneCohesion" parameter="PlaneCohesion"/>
+ <material_property name="PlaneFrictionAngle" parameter="PlaneFrictionAngle"/>
+ <material_property name="PlaneDilatancyAngle" parameter="PlaneDilatancyAngle"/>
+ <material_property name="Cohesion" parameter="Cohesion"/>
+ <material_property name="FrictionAngle" parameter="FrictionAngle"/>
+ <material_property name="DilatancyAngle" parameter="DilatancyAngle"/>
+ <material_property name="TransitionAngle" parameter="TransitionAngle"/>
+ <material_property name="TensionCutOffParameter" parameter="TensionCutOffParameter"/>
+ </material_properties>
+ </constitutive_relation>
+ <solid_density>rho_sr</solid_density>
+ <specific_body_force>0 0</specific_body_force>
+ <process_variables>
+ <process_variable>displacement</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="ElasticStrain" output_name="ElasticStrain"/>
+ <!--secondary_variable type="static" internal_name="EquivalentPlasticStrain" output_name="EquivalentPlasticStrain"/-->
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="SD">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1e-12</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>IterationNumberBasedTimeStepping</type>
+ <t_initial>0.0</t_initial>
+ <t_end>2.0</t_end>
+ <initial_dt>0.01</initial_dt>
+ <minimum_dt>0.001</minimum_dt>
+ <maximum_dt>0.01</maximum_dt>
+ <number_iterations>1 4 10 20 </number_iterations>
+ <multiplier>1.2 1.0 0.9 0.8</multiplier>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>triax_ortho_0_pcs_{:process_id}</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <each_steps>100000</each_steps>
+ </pair>
+ </timesteps>
+ <variables>
+ <variable>displacement</variable>
+ <variable>sigma</variable>
+ </variables>
+ <suffix>_ts_{:timestep}_t_{:time}</suffix>
+ </output>
+ </time_loop>
+ <local_coordinate_system>
+ <basis_vector_0>e0</basis_vector_0>
+ <basis_vector_1>e1</basis_vector_1>
+ </local_coordinate_system>
+ <parameters>
+ <parameter>
+ <name>e0</name>
+ <type>Constant</type>
+ <values>1 0</values>
+ </parameter>
+ <parameter>
+ <name>e1</name>
+ <type>Constant</type>
+ <values>0 1</values>
+ </parameter>
+ <parameter>
+ <name>YoungModulus1</name>
+ <type>Constant</type>
+ <value>29.65e9</value>
+ </parameter>
+ <parameter>
+ <name>YoungModulus2</name>
+ <type>Constant</type>
+ <value>15.2e9</value>
+ </parameter>
+ <parameter>
+ <name>YoungModulus3</name>
+ <type>Constant</type>
+ <value>15.2e9</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio12</name>
+ <type>Constant</type>
+ <value>.2</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio23</name>
+ <type>Constant</type>
+ <value>.22</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio13</name>
+ <type>Constant</type>
+ <value>.2</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus12</name>
+ <type>Constant</type>
+ <value>5.86e9</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus23</name>
+ <type>Constant</type>
+ <value>6.23e9</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus13</name>
+ <type>Constant</type>
+ <value>5.86e9</value>
+ </parameter>
+ <parameter>
+ <name>PlaneCohesion</name>
+ <type>Constant</type>
+ <value>14.e6</value>
+ </parameter>
+ <parameter>
+ <name>PlaneFrictionAngle</name>
+ <type>Constant</type>
+ <value>24</value>
+ </parameter>
+ <parameter>
+ <name>PlaneDilatancyAngle</name>
+ <type>Constant</type>
+ <value>10</value>
+ </parameter>
+ <parameter>
+ <name>Cohesion</name>
+ <type>Constant</type>
+ <value>26.e6</value>
+ </parameter>
+ <parameter>
+ <name>FrictionAngle</name>
+ <type>Constant</type>
+ <value>29</value>
+ </parameter>
+ <parameter>
+ <name>DilatancyAngle</name>
+ <type>Constant</type>
+ <value>10</value>
+ </parameter>
+ <parameter>
+ <name>TransitionAngle</name>
+ <type>Constant</type>
+ <value>29</value>
+ </parameter>
+ <parameter>
+ <name>TensionCutOffParameter</name>
+ <type>Constant</type>
+ <value>1e3</value>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0</values>
+ </parameter>
+ <parameter>
+ <name>zero</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>axial_displacement_top</name>
+ <type>CurveScaled</type>
+ <curve>dis_loading_curve</curve>
+ <parameter>loading_value_top</parameter>
+ </parameter>
+ <parameter>
+ <name>loading_value_top</name>
+ <type>Constant</type>
+ <value>-1e-2</value>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_left</name>
+ <type>Constant</type>
+ <values>1e1</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_left</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_left</parameter>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_right</name>
+ <type>Constant</type>
+ <values>-1e1</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_right</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_right</parameter>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>dis_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 0.0 1.0 </values>
+ </curve>
+ <curve>
+ <name>pre_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 1.0 1.0 </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>square_1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>origin</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>axial_displacement_top</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_right</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_left</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>40</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>SparseLU</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ </eigen>
+ <petsc>
+ <prefix>sd</prefix>
+ <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_1e0_47.prj b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_1e0_47.prj
new file mode 100644
index 0000000000000000000000000000000000000000..00bddb8dca1be42b24d1d8a867e99df61601525c
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_1e0_47.prj
@@ -0,0 +1,324 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>square_1x1_quad_1e0.vtu</mesh>
+ <geometry>square_1x1.gml</geometry>
+ <processes>
+ <process>
+ <name>SD</name>
+ <type>SMALL_DEFORMATION</type>
+ <integration_order>2</integration_order>
+ <constitutive_relation>
+ <type>MFront</type>
+ <behaviour>MohrCoulombAbboSloanAniso</behaviour>
+ <material_properties>
+ <material_property name="YoungModulus1" parameter="YoungModulus1"/>
+ <material_property name="YoungModulus2" parameter="YoungModulus2"/>
+ <material_property name="YoungModulus3" parameter="YoungModulus3"/>
+ <material_property name="PoissonRatio12" parameter="PoissonRatio12"/>
+ <material_property name="PoissonRatio23" parameter="PoissonRatio23"/>
+ <material_property name="PoissonRatio13" parameter="PoissonRatio13"/>
+ <material_property name="ShearModulus12" parameter="ShearModulus12"/>
+ <!--material_property name="ShearModulus23" parameter="ShearModulus23"/>
+ <material_property name="ShearModulus13" parameter="ShearModulus13"/-->
+ <material_property name="Cohesion" parameter="Cohesion"/>
+ <material_property name="FrictionAngle" parameter="FrictionAngle"/>
+ <material_property name="DilatancyAngle" parameter="DilatancyAngle"/>
+ <material_property name="TransitionAngle" parameter="TransitionAngle"/>
+ <material_property name="TensionCutOffParameter" parameter="TensionCutOffParameter"/>
+ <material_property name="NormalFactor" parameter="NormalFactor"/>
+ <material_property name="ShearFactor" parameter="ShearFactor"/>
+ </material_properties>
+ </constitutive_relation>
+ <solid_density>rho_sr</solid_density>
+ <specific_body_force>0 0</specific_body_force>
+ <process_variables>
+ <process_variable>displacement</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="ElasticStrain" output_name="ElasticStrain"/>
+ <!--secondary_variable type="static" internal_name="EquivalentPlasticStrain" output_name="EquivalentPlasticStrain"/-->
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="SD">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1e-14</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>IterationNumberBasedTimeStepping</type>
+ <t_initial>0.0</t_initial>
+ <t_end>2.0</t_end>
+ <initial_dt>0.01</initial_dt>
+ <minimum_dt>0.01</minimum_dt>
+ <maximum_dt>0.01</maximum_dt>
+ <number_iterations>1 4 10 20 </number_iterations>
+ <multiplier>1.2 1.0 0.9 0.8</multiplier>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>triax_aniso_47_pcs_{:process_id}</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <each_steps>100000</each_steps>
+ </pair>
+ </timesteps>
+ <variables>
+ <variable>displacement</variable>
+ <variable>sigma</variable>
+ </variables>
+ <suffix>_ts_{:timestep}_t_{:time}</suffix>
+ </output>
+ </time_loop>
+ <local_coordinate_system>
+ <basis_vector_0>e0</basis_vector_0>
+ <basis_vector_1>e1</basis_vector_1>
+ </local_coordinate_system>
+ <parameters>
+ <parameter>
+ <name>e0</name>
+ <type>Constant</type>
+ <values>0.6862416378687336 0.7273736415730486</values>
+ </parameter>
+ <parameter>
+ <name>e1</name>
+ <type>Constant</type>
+ <values>-0.7273736415730486 0.6862416378687336</values>
+ </parameter>
+ <parameter>
+ <name>YoungModulus1</name>
+ <type>Constant</type>
+ <value>29.65e9</value>
+ </parameter>
+ <parameter>
+ <name>YoungModulus2</name>
+ <type>Constant</type>
+ <value>15.2e9</value>
+ </parameter>
+ <parameter>
+ <name>YoungModulus3</name>
+ <type>Constant</type>
+ <value>15.2e9</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio12</name>
+ <type>Constant</type>
+ <value>.2</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio23</name>
+ <type>Constant</type>
+ <value>.22</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio13</name>
+ <type>Constant</type>
+ <value>.2</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus12</name>
+ <type>Constant</type>
+ <value>5.86e9</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus23</name>
+ <type>Constant</type>
+ <value>6.23e9</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus13</name>
+ <type>Constant</type>
+ <!--value>5.86e9</value-->
+ <value>12.35e9</value>
+ </parameter>
+ <parameter>
+ <name>Cohesion</name>
+ <type>Constant</type>
+ <value>26.e6</value>
+ </parameter>
+ <parameter>
+ <name>FrictionAngle</name>
+ <type>Constant</type>
+ <value>29</value>
+ </parameter>
+ <parameter>
+ <name>DilatancyAngle</name>
+ <type>Constant</type>
+ <value>10</value>
+ </parameter>
+ <parameter>
+ <name>TransitionAngle</name>
+ <type>Constant</type>
+ <value>29</value>
+ </parameter>
+ <parameter>
+ <name>TensionCutOffParameter</name>
+ <type>Constant</type>
+ <value>1e3</value>
+ </parameter>
+ <parameter>
+ <name>ShearFactor</name>
+ <type>Constant</type>
+ <value>1.29</value>
+ </parameter>
+ <parameter>
+ <name>NormalFactor</name>
+ <type>Constant</type>
+ <value>0.92</value>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0</values>
+ </parameter>
+ <parameter>
+ <name>zero</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>axial_displacement_top</name>
+ <type>CurveScaled</type>
+ <curve>dis_loading_curve</curve>
+ <parameter>loading_value_top</parameter>
+ </parameter>
+ <parameter>
+ <name>loading_value_top</name>
+ <type>Constant</type>
+ <value>-1e-2</value>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_left</name>
+ <type>Constant</type>
+ <values>1e1</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_left</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_left</parameter>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_right</name>
+ <type>Constant</type>
+ <values>-1e1</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_right</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_right</parameter>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>dis_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 0.0 1.0 </values>
+ </curve>
+ <curve>
+ <name>pre_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 1.0 1.0 </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>square_1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>origin</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>axial_displacement_top</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_right</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_left</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>40</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>SparseLU</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ </eigen>
+ <petsc>
+ <prefix>sd</prefix>
+ <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+ <test_definition>
+ <vtkdiff>
+ <file>triax_aniso_47_pcs_0_ts_200_t_2.000000.vtu</file>
+ <field>displacement</field>
+ <absolute_tolerance>1e-14</absolute_tolerance>
+ <relative_tolerance>1e-14</relative_tolerance>
+ </vtkdiff>
+ <vtkdiff>
+ <file>triax_aniso_47_pcs_0_ts_200_t_2.000000.vtu</file>
+ <field>sigma</field>
+ <!-- E is 1e10 -> E * u_tol = 1e-4 -->
+ <absolute_tolerance>1e-4</absolute_tolerance>
+ <relative_tolerance>1e-4</relative_tolerance>
+ </vtkdiff>
+ </test_definition>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_47_pcs_0_ts_200_t_2.000000.vtu b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_47_pcs_0_ts_200_t_2.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..2f5f322372ea9a6df6e38942783008b7b8d1e07a
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_47_pcs_0_ts_200_t_2.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ef298fe75f04d4f006a5c2dc594848eb8aa6bb78e4389b112ba548aa18d91dd6
+size 5324
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_original.prj b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_original.prj
new file mode 100644
index 0000000000000000000000000000000000000000..e39250e98b81e2373af42a741e05227215ed0471
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_aniso_original.prj
@@ -0,0 +1,309 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>square_1x1_quad_1e0.vtu</mesh>
+ <geometry>square_1x1.gml</geometry>
+ <processes>
+ <process>
+ <name>SD</name>
+ <type>SMALL_DEFORMATION</type>
+ <integration_order>2</integration_order>
+ <constitutive_relation>
+ <type>MFront</type>
+ <behaviour>MohrCoulombAbboSloanAniso</behaviour>
+ <material_properties>
+ <material_property name="YoungModulus1" parameter="YoungModulus1"/>
+ <material_property name="YoungModulus2" parameter="YoungModulus2"/>
+ <material_property name="YoungModulus3" parameter="YoungModulus3"/>
+ <material_property name="PoissonRatio12" parameter="PoissonRatio12"/>
+ <material_property name="PoissonRatio23" parameter="PoissonRatio23"/>
+ <material_property name="PoissonRatio13" parameter="PoissonRatio13"/>
+ <material_property name="ShearModulus12" parameter="ShearModulus12"/>
+ <!--material_property name="ShearModulus23" parameter="ShearModulus23"/>
+ <material_property name="ShearModulus13" parameter="ShearModulus13"/-->
+ <material_property name="Cohesion" parameter="Cohesion"/>
+ <material_property name="FrictionAngle" parameter="FrictionAngle"/>
+ <material_property name="DilatancyAngle" parameter="DilatancyAngle"/>
+ <material_property name="TransitionAngle" parameter="TransitionAngle"/>
+ <material_property name="TensionCutOffParameter" parameter="TensionCutOffParameter"/>
+ <material_property name="NormalFactor" parameter="NormalFactor"/>
+ <material_property name="ShearFactor" parameter="ShearFactor"/>
+ </material_properties>
+ </constitutive_relation>
+ <solid_density>rho_sr</solid_density>
+ <specific_body_force>0 0</specific_body_force>
+ <process_variables>
+ <process_variable>displacement</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="ElasticStrain" output_name="ElasticStrain"/>
+ <!--secondary_variable type="static" internal_name="EquivalentPlasticStrain" output_name="EquivalentPlasticStrain"/-->
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="SD">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1e-12</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>IterationNumberBasedTimeStepping</type>
+ <t_initial>0.0</t_initial>
+ <t_end>2.0</t_end>
+ <initial_dt>0.01</initial_dt>
+ <minimum_dt>0.001</minimum_dt>
+ <maximum_dt>0.01</maximum_dt>
+ <number_iterations>1 4 10 20 </number_iterations>
+ <multiplier>1.2 1.0 0.9 0.8</multiplier>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>triax_aniso_0_pcs_{:process_id}</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <each_steps>100000</each_steps>
+ </pair>
+ </timesteps>
+ <variables>
+ <variable>displacement</variable>
+ <variable>sigma</variable>
+ </variables>
+ <suffix>_ts_{:timestep}_t_{:time}</suffix>
+ </output>
+ </time_loop>
+ <local_coordinate_system>
+ <basis_vector_0>e0</basis_vector_0>
+ <basis_vector_1>e1</basis_vector_1>
+ </local_coordinate_system>
+ <parameters>
+ <parameter>
+ <name>e0</name>
+ <type>Constant</type>
+ <values>1 0</values>
+ </parameter>
+ <parameter>
+ <name>e1</name>
+ <type>Constant</type>
+ <values>0 1</values>
+ </parameter>
+ <parameter>
+ <name>YoungModulus1</name>
+ <type>Constant</type>
+ <value>29.65e9</value>
+ </parameter>
+ <parameter>
+ <name>YoungModulus2</name>
+ <type>Constant</type>
+ <value>15.2e9</value>
+ </parameter>
+ <parameter>
+ <name>YoungModulus3</name>
+ <type>Constant</type>
+ <value>15.2e9</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio12</name>
+ <type>Constant</type>
+ <value>.2</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio23</name>
+ <type>Constant</type>
+ <value>.22</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio13</name>
+ <type>Constant</type>
+ <value>.2</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus12</name>
+ <type>Constant</type>
+ <value>5.86e9</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus23</name>
+ <type>Constant</type>
+ <value>6.23e9</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus13</name>
+ <type>Constant</type>
+ <!--value>5.86e9</value-->
+ <value>12.35e9</value>
+ </parameter>
+ <parameter>
+ <name>Cohesion</name>
+ <type>Constant</type>
+ <value>26.e6</value>
+ </parameter>
+ <parameter>
+ <name>FrictionAngle</name>
+ <type>Constant</type>
+ <value>29</value>
+ </parameter>
+ <parameter>
+ <name>DilatancyAngle</name>
+ <type>Constant</type>
+ <value>10</value>
+ </parameter>
+ <parameter>
+ <name>TransitionAngle</name>
+ <type>Constant</type>
+ <value>29</value>
+ </parameter>
+ <parameter>
+ <name>TensionCutOffParameter</name>
+ <type>Constant</type>
+ <value>1e3</value>
+ </parameter>
+ <parameter>
+ <name>ShearFactor</name>
+ <type>Constant</type>
+ <value>1.29</value>
+ </parameter>
+ <parameter>
+ <name>NormalFactor</name>
+ <type>Constant</type>
+ <value>0.92</value>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0</values>
+ </parameter>
+ <parameter>
+ <name>zero</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>axial_displacement_top</name>
+ <type>CurveScaled</type>
+ <curve>dis_loading_curve</curve>
+ <parameter>loading_value_top</parameter>
+ </parameter>
+ <parameter>
+ <name>loading_value_top</name>
+ <type>Constant</type>
+ <value>-1e-2</value>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_left</name>
+ <type>Constant</type>
+ <values>1e1</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_left</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_left</parameter>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_right</name>
+ <type>Constant</type>
+ <values>-1e1</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_right</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_right</parameter>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>dis_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 0.0 1.0 </values>
+ </curve>
+ <curve>
+ <name>pre_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 1.0 1.0 </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>square_1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>origin</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>axial_displacement_top</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_right</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_left</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>40</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>SparseLU</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ </eigen>
+ <petsc>
+ <prefix>sd</prefix>
+ <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_original.prj b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_original.prj
new file mode 100644
index 0000000000000000000000000000000000000000..b57771f060496adf7a8d23f7f14e59604f1c596d
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_original.prj
@@ -0,0 +1,272 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>square_1x1_quad_1e0.vtu</mesh>
+ <geometry>square_1x1.gml</geometry>
+ <processes>
+ <process>
+ <name>SD</name>
+ <type>SMALL_DEFORMATION</type>
+ <integration_order>2</integration_order>
+ <constitutive_relation>
+ <type>MFront</type>
+ <behaviour>MohrCoulombAbboSloanUBI</behaviour>
+ <material_properties>
+ <material_property name="YoungModulus" parameter="YoungModulus"/>
+ <material_property name="PoissonRatio" parameter="PoissonRatio"/>
+ <material_property name="PlaneCohesion" parameter="PlaneCohesion"/>
+ <material_property name="PlaneFrictionAngle" parameter="PlaneFrictionAngle"/>
+ <material_property name="PlaneDilatancyAngle" parameter="PlaneDilatancyAngle"/>
+ <material_property name="Cohesion" parameter="Cohesion"/>
+ <material_property name="FrictionAngle" parameter="FrictionAngle"/>
+ <material_property name="DilatancyAngle" parameter="DilatancyAngle"/>
+ <material_property name="TransitionAngle" parameter="TransitionAngle"/>
+ <material_property name="TensionCutOffParameter" parameter="TensionCutOffParameter"/>
+ </material_properties>
+ </constitutive_relation>
+ <solid_density>rho_sr</solid_density>
+ <specific_body_force>0 0</specific_body_force>
+ <process_variables>
+ <process_variable>displacement</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="ElasticStrain" output_name="ElasticStrain"/>
+ <!--secondary_variable type="static" internal_name="EquivalentPlasticStrain" output_name="EquivalentPlasticStrain"/-->
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="SD">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1e-12</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>IterationNumberBasedTimeStepping</type>
+ <t_initial>0.0</t_initial>
+ <t_end>2.0</t_end>
+ <initial_dt>0.01</initial_dt>
+ <minimum_dt>0.001</minimum_dt>
+ <maximum_dt>0.01</maximum_dt>
+ <number_iterations>1 4 10 20 </number_iterations>
+ <multiplier>1.2 1.0 0.9 0.8</multiplier>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>triax_0_pcs_{:process_id}</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <each_steps>10000</each_steps>
+ </pair>
+ </timesteps>
+ <variables>
+ <variable>displacement</variable>
+ <variable>sigma</variable>
+ </variables>
+ <suffix>_ts_{:timestep}_t_{:time}</suffix>
+ </output>
+ </time_loop>
+ <local_coordinate_system>
+ <basis_vector_0>e0</basis_vector_0>
+ <basis_vector_1>e1</basis_vector_1>
+ </local_coordinate_system>
+ <parameters>
+ <parameter>
+ <name>e0</name>
+ <type>Constant</type>
+ <values>1 0</values>
+ </parameter>
+ <parameter>
+ <name>e1</name>
+ <type>Constant</type>
+ <values>0 1</values>
+ </parameter>
+ <parameter>
+ <name>YoungModulus</name>
+ <type>Constant</type>
+ <value>170e6</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio</name>
+ <type>Constant</type>
+ <value>.22</value>
+ </parameter>
+ <parameter>
+ <name>PlaneCohesion</name>
+ <type>Constant</type>
+ <value>1e3</value>
+ </parameter>
+ <parameter>
+ <name>PlaneFrictionAngle</name>
+ <type>Constant</type>
+ <value>30</value>
+ </parameter>
+ <parameter>
+ <name>PlaneDilatancyAngle</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>Cohesion</name>
+ <type>Constant</type>
+ <value>2.e3</value>
+ </parameter>
+ <parameter>
+ <name>FrictionAngle</name>
+ <type>Constant</type>
+ <value>40</value>
+ </parameter>
+ <parameter>
+ <name>DilatancyAngle</name>
+ <type>Constant</type>
+ <value>10</value>
+ </parameter>
+ <parameter>
+ <name>TransitionAngle</name>
+ <type>Constant</type>
+ <value>29</value>
+ </parameter>
+ <parameter>
+ <name>TensionCutOffParameter</name>
+ <type>Constant</type>
+ <value>1.0</value>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0</values>
+ </parameter>
+ <parameter>
+ <name>zero</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>axial_displacement_top</name>
+ <type>CurveScaled</type>
+ <curve>dis_loading_curve</curve>
+ <parameter>loading_value_top</parameter>
+ </parameter>
+ <parameter>
+ <name>loading_value_top</name>
+ <type>Constant</type>
+ <value>-1.e-4</value>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_left</name>
+ <type>Constant</type>
+ <values>1e0</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_left</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_left</parameter>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_right</name>
+ <type>Constant</type>
+ <values>-1e0</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_right</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_right</parameter>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>dis_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 0.0 1.0 </values>
+ </curve>
+ <curve>
+ <name>pre_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 1.0 1.0 </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>square_1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>origin</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>axial_displacement_top</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_right</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_left</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>40</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>SparseLU</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ </eigen>
+ <petsc>
+ <prefix>sd</prefix>
+ <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ortho_1e0_47.prj b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ortho_1e0_47.prj
new file mode 100644
index 0000000000000000000000000000000000000000..c7c8d7bb19ad80a6310d8b5e017f089212640a63
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ortho_1e0_47.prj
@@ -0,0 +1,329 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>square_1x1_quad_1e0.vtu</mesh>
+ <geometry>square_1x1.gml</geometry>
+ <processes>
+ <process>
+ <name>SD</name>
+ <type>SMALL_DEFORMATION</type>
+ <integration_order>2</integration_order>
+ <constitutive_relation>
+ <type>MFront</type>
+ <behaviour>MohrCoulombAbboSloanUBIOrtho</behaviour>
+ <material_properties>
+ <material_property name="YoungModulus1" parameter="YoungModulus1"/>
+ <material_property name="YoungModulus2" parameter="YoungModulus2"/>
+ <material_property name="YoungModulus3" parameter="YoungModulus3"/>
+ <material_property name="PoissonRatio12" parameter="PoissonRatio12"/>
+ <material_property name="PoissonRatio23" parameter="PoissonRatio23"/>
+ <material_property name="PoissonRatio13" parameter="PoissonRatio13"/>
+ <material_property name="ShearModulus12" parameter="ShearModulus12"/>
+ <!--material_property name="ShearModulus23" parameter="ShearModulus23"/>
+ <material_property name="ShearModulus13" parameter="ShearModulus13"/-->
+ <material_property name="PlaneCohesion" parameter="PlaneCohesion"/>
+ <material_property name="PlaneFrictionAngle" parameter="PlaneFrictionAngle"/>
+ <material_property name="PlaneDilatancyAngle" parameter="PlaneDilatancyAngle"/>
+ <material_property name="Cohesion" parameter="Cohesion"/>
+ <material_property name="FrictionAngle" parameter="FrictionAngle"/>
+ <material_property name="DilatancyAngle" parameter="DilatancyAngle"/>
+ <material_property name="TransitionAngle" parameter="TransitionAngle"/>
+ <material_property name="TensionCutOffParameter" parameter="TensionCutOffParameter"/>
+ </material_properties>
+ </constitutive_relation>
+ <solid_density>rho_sr</solid_density>
+ <specific_body_force>0 0</specific_body_force>
+ <process_variables>
+ <process_variable>displacement</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="ElasticStrain" output_name="ElasticStrain"/>
+ <!--secondary_variable type="static" internal_name="EquivalentPlasticStrain" output_name="EquivalentPlasticStrain"/-->
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="SD">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1e-14</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>IterationNumberBasedTimeStepping</type>
+ <t_initial>0.0</t_initial>
+ <t_end>2.0</t_end>
+ <initial_dt>0.01</initial_dt>
+ <minimum_dt>0.01</minimum_dt>
+ <maximum_dt>0.01</maximum_dt>
+ <number_iterations>1 4 10 20 </number_iterations>
+ <multiplier>1.2 1.0 0.9 0.8</multiplier>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>triax_ortho_47_pcs_{:process_id}</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <each_steps>100000</each_steps>
+ </pair>
+ </timesteps>
+ <variables>
+ <variable>displacement</variable>
+ <variable>sigma</variable>
+ </variables>
+ <suffix>_ts_{:timestep}_t_{:time}</suffix>
+ </output>
+ </time_loop>
+ <local_coordinate_system>
+ <basis_vector_0>e0</basis_vector_0>
+ <basis_vector_1>e1</basis_vector_1>
+ </local_coordinate_system>
+ <parameters>
+ <parameter>
+ <name>e0</name>
+ <type>Constant</type>
+ <values>0.6862416378687336 0.7273736415730486</values>
+ </parameter>
+ <parameter>
+ <name>e1</name>
+ <type>Constant</type>
+ <values>-0.7273736415730486 0.6862416378687336</values>
+ </parameter>
+ <parameter>
+ <name>YoungModulus1</name>
+ <type>Constant</type>
+ <value>29.65e9</value>
+ </parameter>
+ <parameter>
+ <name>YoungModulus2</name>
+ <type>Constant</type>
+ <value>15.2e9</value>
+ </parameter>
+ <parameter>
+ <name>YoungModulus3</name>
+ <type>Constant</type>
+ <value>15.2e9</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio12</name>
+ <type>Constant</type>
+ <value>.2</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio23</name>
+ <type>Constant</type>
+ <value>.22</value>
+ </parameter>
+ <parameter>
+ <name>PoissonRatio13</name>
+ <type>Constant</type>
+ <value>.2</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus12</name>
+ <type>Constant</type>
+ <value>5.86e9</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus23</name>
+ <type>Constant</type>
+ <value>6.23e9</value>
+ </parameter>
+ <parameter>
+ <name>ShearModulus13</name>
+ <type>Constant</type>
+ <value>5.86e9</value>
+ </parameter>
+ <parameter>
+ <name>PlaneCohesion</name>
+ <type>Constant</type>
+ <value>14.e6</value>
+ </parameter>
+ <parameter>
+ <name>PlaneFrictionAngle</name>
+ <type>Constant</type>
+ <value>24</value>
+ </parameter>
+ <parameter>
+ <name>PlaneDilatancyAngle</name>
+ <type>Constant</type>
+ <value>10</value>
+ </parameter>
+ <parameter>
+ <name>Cohesion</name>
+ <type>Constant</type>
+ <value>26.e6</value>
+ </parameter>
+ <parameter>
+ <name>FrictionAngle</name>
+ <type>Constant</type>
+ <value>29</value>
+ </parameter>
+ <parameter>
+ <name>DilatancyAngle</name>
+ <type>Constant</type>
+ <value>10</value>
+ </parameter>
+ <parameter>
+ <name>TransitionAngle</name>
+ <type>Constant</type>
+ <value>29</value>
+ </parameter>
+ <parameter>
+ <name>TensionCutOffParameter</name>
+ <type>Constant</type>
+ <value>1e3</value>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0</values>
+ </parameter>
+ <parameter>
+ <name>zero</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>axial_displacement_top</name>
+ <type>CurveScaled</type>
+ <curve>dis_loading_curve</curve>
+ <parameter>loading_value_top</parameter>
+ </parameter>
+ <parameter>
+ <name>loading_value_top</name>
+ <type>Constant</type>
+ <value>-1e-2</value>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_left</name>
+ <type>Constant</type>
+ <values>1e1</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_left</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_left</parameter>
+ </parameter>
+ <parameter>
+ <name>neumann_force_value_right</name>
+ <type>Constant</type>
+ <values>-1e1</values>
+ </parameter>
+ <parameter>
+ <name>neumann_force_right</name>
+ <type>CurveScaled</type>
+ <curve>pre_loading_curve</curve>
+ <parameter>neumann_force_value_right</parameter>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>dis_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 0.0 1.0 </values>
+ </curve>
+ <curve>
+ <name>pre_loading_curve</name>
+ <coords>0.0 1.0 2.0 </coords>
+ <values>0.0 1.0 1.0 </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>square_1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>origin</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>zero</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>axial_displacement_top</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_right</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Neumann</type>
+ <component>0</component>
+ <parameter>neumann_force_left</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>40</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>SparseLU</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ </eigen>
+ <petsc>
+ <prefix>sd</prefix>
+ <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+ <test_definition>
+ <vtkdiff>
+ <file>triax_ortho_47_pcs_0_ts_200_t_2.000000.vtu</file>
+ <field>displacement</field>
+ <absolute_tolerance>1e-14</absolute_tolerance>
+ <relative_tolerance>1e-14</relative_tolerance>
+ </vtkdiff>
+ <vtkdiff>
+ <file>triax_ortho_47_pcs_0_ts_200_t_2.000000.vtu</file>
+ <field>sigma</field>
+ <!-- E is 1e10 -> E * u_tol = 1e-4 -->
+ <absolute_tolerance>1e-4</absolute_tolerance>
+ <relative_tolerance>1e-4</relative_tolerance>
+ </vtkdiff>
+ </test_definition>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ortho_47_pcs_0_ts_200_t_2.000000.vtu b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ortho_47_pcs_0_ts_200_t_2.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..85cde96dc1db1514091473c8bbb418abfd2a4875
--- /dev/null
+++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/triax_ortho_47_pcs_0_ts_200_t_2.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9ffbd870c833690c90a5667ce787ac55cab3e8d34acd8c65bcc9a57692807f31
+size 5328