Bug fix in J2 plasticity routine

FANS_Dashboard/fans_dashboard/plotting/plotting.py

@@ -149,8 +149,8 @@ def plot_subplots(
 
     # Update layout with the overall plot title and styling
     fig.update_layout(
-        height=1000,
-        width=1600,
+        height=500,
+        width=800,
         title_text=title,
         title_font=dict(size=fontsize),
         showlegend=False,  # Legends removed

FANS_Dashboard/pyproject.toml

@@ -13,6 +13,7 @@ dependencies = [
   "lxml",
   "scipy",
   "meshio",
+  "nbformat",
 ]
 
 [tool.hatch.build.targets.wheel]

include/material_models/J2Plasticity.h

@@ -22,17 +22,12 @@ class J2Plasticity : public MechModel {
         }
         n_mat = bulk_modulus.size();
 
-        Matrix<double, 6, 6> *Ce      = new Matrix<double, 6, 6>[n_mat];
-        Matrix<double, 6, 6>  topLeft = Matrix<double, 6, 6>::Zero();
-        topLeft.topLeftCorner(3, 3).setConstant(1);
-
-        kapparef_mat = Matrix<double, n_str, n_str>::Zero();
-        for (int i = 0; i < n_mat; ++i) {
-            Ce[i] = 3 * bulk_modulus[i] * topLeft +
-                    2 * shear_modulus[i] * (-1.0 / 3.0 * topLeft + Matrix<double, 6, 6>::Identity());
-            kapparef_mat += Ce[i];
-        }
-        kapparef_mat /= n_mat;
+        const double Kbar      = std::accumulate(bulk_modulus.begin(), bulk_modulus.end(), 0.0) / static_cast<double>(n_mat);
+        const double Gbar      = std::accumulate(shear_modulus.begin(), shear_modulus.end(), 0.0) / static_cast<double>(n_mat);
+        const double lambdabar = Kbar - 2.0 * Gbar / 3.0;
+        kapparef_mat.setZero();
+        kapparef_mat.topLeftCorner(3, 3).setConstant(lambdabar);
+        kapparef_mat.diagonal().array() += 2.0 * Gbar;
 
         // Allocate the member matrices/vectors for performance optimization
         sqrt_two_over_three = sqrt(2.0 / 3.0);
@@ -81,7 +76,7 @@ class J2Plasticity : public MechModel {
         treps       = eps_elastic.head<3>().sum();
 
         // Compute trial stress
-        sigma_trial_n1.head<3>().setConstant(bulk_modulus[mat_index] * treps);
+        sigma_trial_n1.head<3>().setConstant((bulk_modulus[mat_index] - 2.0 * shear_modulus[mat_index] / 3.0) * treps);
         sigma_trial_n1.head<3>() += 2 * shear_modulus[mat_index] * eps_elastic.head<3>();
         sigma_trial_n1.tail<3>() = 2 * shear_modulus[mat_index] * eps_elastic.tail<3>();
 
@@ -90,16 +85,19 @@ class J2Plasticity : public MechModel {
         dev.head<3>().array() -= sigma_trial_n1.head<3>().mean();
 
         // Compute trial q and q_bar
-        q_trial_n1              = compute_q_trial(psi_t[element_idx](i / n_str), mat_index);
-        qbar_trial_n1.head<3>() = -H[mat_index] * (2.0 / 3.0) * psi_bar_t[element_idx].col(i / n_str).head<3>();
-        qbar_trial_n1.tail<3>().setZero(); // Lower part is zero
+        q_trial_n1    = compute_q_trial(psi_t[element_idx](i / n_str), mat_index);
+        qbar_trial_n1 = -(2.0 / 3.0) * H[mat_index] * psi_bar_t[element_idx].col(i / n_str);
 
         // Calculate the trial yield function
         dev_minus_qbar      = dev - qbar_trial_n1;
         norm_dev_minus_qbar = dev_minus_qbar.norm();
 
         // Avoid division by zero
-        n = (norm_dev_minus_qbar < 1e-12) ? n.setZero() : dev_minus_qbar / norm_dev_minus_qbar;
+        if (norm_dev_minus_qbar < 1e-12) {
+            n.setZero();
+        } else {
+            n = dev_minus_qbar / norm_dev_minus_qbar;
+        }
 
         f_trial = norm_dev_minus_qbar - sqrt_two_over_three * (yield_stress[mat_index] - q_trial_n1);
 
@@ -191,7 +189,7 @@ class J2ViscoPlastic_NonLinearIsotropicHardening : public J2Plasticity {
         denominator.resize(n_mat);
         sigma_diff.resize(n_mat);
         for (size_t i = 0; i < n_mat; ++i) {
-            denominator[i] = 2 * shear_modulus[i] + H[i] * (2.0 / 3.0) + eta[i] / dt;
+            denominator[i] = 2 * shear_modulus[i] + (2.0 / 3.0) * (K[i] + H[i]) + eta[i] / dt;
             sigma_diff[i]  = sqrt_two_over_three * (sigma_inf[i] - yield_stress[i]);
         }
     }
@@ -239,7 +237,7 @@ class J2ViscoPlastic_NonLinearIsotropicHardening : public J2Plasticity {
     vector<double> sigma_diff;  // sqrt(2/3) * (sigma_inf - yield_stress)
 };
 
-void J2Plasticity::postprocess(Solver<3> &solver, Reader &reader, const char *resultsFileName, int load_idx, int time_idx)
+inline void J2Plasticity::postprocess(Solver<3> &solver, Reader &reader, const char *resultsFileName, int load_idx, int time_idx)
 {
     int      n_str                             = 6; // The plastic strain and stress vectors have 6 components each
     VectorXd mean_plastic_strain               = VectorXd::Zero(solver.local_n0 * solver.n_y * solver.n_z * n_str);

include/material_models/PseudoPlastic.h

@@ -32,18 +32,12 @@ class PseudoPlastic : public MechModel {
         }
         n_mat = bulk_modulus.size();
 
-        // Initialize stiffness matrix (assuming for two materials, otherwise needs extension)
-        Matrix<double, 6, 6> *Ce      = new Matrix<double, 6, 6>[n_mat];
-        Matrix<double, 6, 6>  topLeft = Matrix<double, 6, 6>::Zero();
-        topLeft.topLeftCorner(3, 3).setConstant(1);
-
-        kapparef_mat = Matrix<double, n_str, n_str>::Zero();
-        for (int i = 0; i < n_mat; ++i) {
-            Ce[i] = 3 * bulk_modulus[i] * topLeft +
-                    2 * shear_modulus[i] * (-1.0 / 3.0 * topLeft + Matrix<double, 6, 6>::Identity());
-            kapparef_mat += Ce[i];
-        }
-        kapparef_mat /= n_mat;
+        const double Kbar      = std::accumulate(bulk_modulus.begin(), bulk_modulus.end(), 0.0) / static_cast<double>(n_mat);
+        const double Gbar      = std::accumulate(shear_modulus.begin(), shear_modulus.end(), 0.0) / static_cast<double>(n_mat);
+        const double lambdabar = Kbar - 2.0 * Gbar / 3.0;
+        kapparef_mat.setZero();
+        kapparef_mat.topLeftCorner(3, 3).setConstant(lambdabar);
+        kapparef_mat.diagonal().array() += 2.0 * Gbar;
     }
 
     void initializeInternalVariables(ptrdiff_t num_elements, int num_gauss_points) override