[Mat/Fracture] support 3d in MC model

MaterialLib/FractureModels/MohrCoulomb.cpp

@@ -56,24 +56,23 @@ void MohrCoulomb<DisplacementDim>::computeConstitutiveRelation(
         typename FractureModelBase<DisplacementDim>::MaterialStateVariables&
         material_state_variables)
 {
-    if (DisplacementDim == 3)
-    {
-        OGS_FATAL("MohrCoulomb fracture model does not support 3D case.");
-        return;
-    }
     material_state_variables.reset();
 
     MaterialPropertyValues const mat(_mp, t, x);
     Eigen::VectorXd const dw = w - w_prev;
 
-    Eigen::MatrixXd Ke = Eigen::MatrixXd::Zero(2, 2);
-    Ke(0,0) = mat.Ks;
-    Ke(1,1) = mat.Kn;
+    const int index_ns = DisplacementDim - 1;
+    Eigen::MatrixXd Ke = Eigen::MatrixXd::Zero(DisplacementDim,DisplacementDim);
+    for (int i=0; i<index_ns; i++)
+        Ke(i,i) = mat.Ks;
+    Ke(index_ns, index_ns) = mat.Kn;
 
     sigma.noalias() = sigma_prev + Ke * dw;
+    double const tau = (DisplacementDim==2) ? sigma[0] : std::sqrt(sigma[0]*sigma[0] + sigma[1]*sigma[1]);
+    double const sigma_n = sigma[index_ns];
 
     // if opening
-    if (sigma[1] > 0)
+    if (sigma_n > 0)
     {
         Kep.setZero();
         sigma.setZero();
@@ -82,7 +81,8 @@ void MohrCoulomb<DisplacementDim>::computeConstitutiveRelation(
     }
 
     // check shear yield function (Fs)
-    double const Fs = std::abs(sigma[0]) + sigma[1] * std::tan(mat.phi) - mat.c;
+    double const Fs = std::abs(tau) + sigma_n * std::tan(mat.phi) - mat.c;
+
     material_state_variables.setShearYieldFunctionValue(Fs);
     if (Fs < .0)
     {
@@ -90,14 +90,25 @@ void MohrCoulomb<DisplacementDim>::computeConstitutiveRelation(
         return;
     }
 
-    Eigen::VectorXd dFs_dS(2);
-    dFs_dS[0] = boost::math::sign(sigma[0]);
-    dFs_dS[1] = std::tan(mat.phi);
-
+    Eigen::VectorXd dFs_dS(DisplacementDim);
     // plastic potential function: Qs = |tau| + Sn * tan da
-    Eigen::VectorXd dQs_dS(2);
-    dQs_dS[0] = boost::math::sign(sigma[0]);
-    dQs_dS[1] = std::tan(mat.psi);
+    Eigen::VectorXd dQs_dS(DisplacementDim);
+
+    if (DisplacementDim == 2)
+    {
+        dFs_dS[0] = boost::math::sign(tau);
+        dQs_dS[0] = boost::math::sign(tau);
+    }
+    else
+    {
+        for (int i=0; i<index_ns-1; i++)
+            dFs_dS[i] = boost::math::sign(tau) / tau * sigma[i];
+
+        for (int i=0; i<index_ns-1; i++)
+            dQs_dS[i] = boost::math::sign(tau) / tau * sigma[i];
+    }
+    dFs_dS[index_ns] = std::tan(mat.phi);
+    dQs_dS[index_ns] = std::tan(mat.psi);
 
     // plastic multiplier
     Eigen::RowVectorXd const A = dFs_dS.transpose() * Ke / (dFs_dS.transpose() * Ke * dQs_dS);