Add finite strain J2 plasticity model

.gitignore

@@ -209,6 +209,7 @@ test/input_files/**/*.json
 !test_J2Plasticity.json
 !test_MixedBCs.json
 !test_CompressibleNeoHookean.json
+!test_FiniteStrainJ2Plasticity.json
 !test_MixedBCs_LargeStrain.json
 
 # pixi environments

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## latest
 
+- Add rate-independent finite-strain J2 plasticity with linear isotropic hardening [#149](https://github.com/DataAnalyticsEngineering/FANS/pull/149)
 - Added large strain support for pyFANS [#141](https://github.com/DataAnalyticsEngineering/FANS/pull/141)
 - Bugfix: fixes GBDiffusion post MaterialManager update and adds a test case [#147](https://github.com/DataAnalyticsEngineering/FANS/pull/147)
 - Bugfix: mem-leak while reading microstructure and added wider CMake support [#140](https://github.com/DataAnalyticsEngineering/FANS/pull/140)

CMakeLists.txt

@@ -151,6 +151,7 @@ target_include_directories(FANS_FANS PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SO
 set_property(TARGET FANS_FANS PROPERTY PUBLIC_HEADER
         include/general.h
         include/matmodel.h
+        include/LargeStrainMechModel.h
         include/MaterialManager.h
         include/reader.h
         include/solverCG.h
@@ -159,16 +160,17 @@ set_property(TARGET FANS_FANS PROPERTY PUBLIC_HEADER
         include/setup.h
         include/mixedBCs.h
 
-        include/material_models/LinearThermal.h
-        include/material_models/GBDiffusion.h
+        include/material_models/thermal/LinearThermal.h
+        include/material_models/thermal/GBDiffusion.h
 
-        include/material_models/LinearElastic.h
-        include/material_models/PseudoPlastic.h
-        include/material_models/J2Plasticity.h
-        include/material_models/J2PlasticityNew.h
+        include/material_models/small_strain/LinearElastic.h
+        include/material_models/small_strain/PseudoPlastic.h
+        include/material_models/small_strain/J2Plasticity.h
+        include/material_models/small_strain/J2PlasticityNew.h
 
-        include/material_models/SaintVenantKirchhoff.h
-        include/material_models/CompressibleNeoHookean.h
+        include/material_models/large_strain/SaintVenantKirchhoff.h
+        include/material_models/large_strain/CompressibleNeoHookean.h
+        include/material_models/large_strain/FiniteStrainJ2Plasticity.h
 )
 
 # version.h is generated and added to the build include directory

README.md

@@ -279,6 +279,7 @@ FANS requires a JSON input file specifying the problem parameters. Example input
     - `J2ViscoPlastic_LinearIsotropicHardening` / `J2ViscoPlastic_NonLinearIsotropicHardening` for rate-independent / dependent J2 plasticity model with kinematic and linear/nonlinear isotropic hardening.
     - `SaintVenantKirchhoff` for the hyperelastic Saint Venant-Kirchhoff material model.
     - `CompressibleNeoHookean` for the compressible Neo-Hookean material model.
+    - `FiniteStrainJ2Plasticity` for rate-independent finite-strain J2 plasticity with linear isotropic hardening.
 
   - `material_properties`: Material parameters specific to the chosen model. Properties are defined as arrays, where each element corresponds to one of the phases listed in the `phases` array.
 

include/LargeStrainMechModel.h

@@ -6,7 +6,7 @@
 /**
  * @brief Base class for large strain mechanical material models
  * This class provides the kinematic framework for finite deformation mechanics.
- * Derived classes only need to implement the constitutive response (S from E).
+ * Derived classes may implement the constitutive response through either compute_S() or get_sigma().
  */
 class LargeStrainMechModel : public Matmodel<3, 9> {
   public:
@@ -74,18 +74,20 @@ class LargeStrainMechModel : public Matmodel<3, 9> {
     // ============================================================================
 
     /**
-     * @brief Compute 2nd Piola-Kirchhoff stress and material tangent C = dS/dE from deformation gradient
-     *
-     * This is the core constitutive function that derived classes must implement.
-     * It should compute S based on the material's constitutive law.
-     * @param F Deformation gradient (3×3)
-     * @param mat_index Material index
-     * @param element_idx Element index
-     * @return 2nd Piola-Kirchhoff stress S (symmetric 3×3)
-     * @return the 6×6 material tangent C = dS/dE in Mandel notation: Ordering: (11, 22, 33, sqrt(2)·12, sqrt(2)·13, sqrt(2)·23)
+     * @brief Convenience constitutive interface for models formulated in terms of S
+     * @param F Deformation gradient
+     * @param mat_index Material-property index local to this material model
+     * @param element_idx Local element index
+     * @param i Offset of this Gauss point in the flattened strain vector
      */
-    virtual Matrix3d             compute_S(const Matrix3d &F, int mat_index, ptrdiff_t element_idx) = 0;
-    virtual Matrix<double, 6, 6> compute_material_tangent(const Matrix3d &F, int mat_index)         = 0;
+    virtual Matrix3d compute_S(const Matrix3d &F, int mat_index, ptrdiff_t element_idx, int i)
+    {
+        throw std::logic_error("Large-strain material model must override compute_S() or get_sigma().");
+    }
+    virtual Matrix<double, 6, 6> compute_material_tangent(const Matrix3d &F, int mat_index, ptrdiff_t element_idx, int i)
+    {
+        throw std::logic_error("Algorithmic material tangent is not implemented for this large-strain model.");
+    }
 
     /**
      * @brief Convert material tangent C (dS/dE) to spatial tangent A (dP/dF)
@@ -190,19 +192,19 @@ class LargeStrainMechModel : public Matmodel<3, 9> {
     Matrix<double, 9, 24> Compute_B(const double x, const double y, const double z) override;
 
     /**
-     * @brief Compute stress at Gauss point (implements base class interface)
+     * @brief Default stress update for models formulated in terms of S
      * This function:
      * 1. Extracts F from eps
      * 2. Calls compute_S(F) - material model implements this
      * 3. Computes P = F S
      * 4. Stores P in sigma
      */
-    void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override final
+    void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override
     {
-        Matrix3d F = extract_F(i);                         // Extract deformation gradient
-        Matrix3d S = compute_S(F, mat_index, element_idx); // Material model computes 2nd PK stress from F
-        Matrix3d P = push_forward(F, S);                   // Push forward to 1st PK stress
-        store_P(i, P);                                     // Store in sigma array
+        Matrix3d F = extract_F(i);                            // Extract deformation gradient
+        Matrix3d S = compute_S(F, mat_index, element_idx, i); // Material model computes 2nd PK stress from F
+        Matrix3d P = push_forward(F, S);                      // Push forward to 1st PK stress
+        store_P(i, P);                                        // Store in sigma array
     }
 };
 

include/material_models/CompressibleNeoHookean.h → include/material_models/large_strain/CompressibleNeoHookean.h

@@ -32,7 +32,7 @@ class CompressibleNeoHookean : public LargeStrainMechModel {
         }
     }
 
-    Matrix3d compute_S(const Matrix3d &F, int mat_index, ptrdiff_t element_idx) override
+    Matrix3d compute_S(const Matrix3d &F, int mat_index, ptrdiff_t element_idx, int i) override
     {
         Matrix3d C = compute_C(F);
         double   J = F.determinant();
@@ -47,7 +47,7 @@ class CompressibleNeoHookean : public LargeStrainMechModel {
         return lambda[mat_index] * logJ * C_inv + mu[mat_index] * (Matrix3d::Identity() - C_inv);
     }
 
-    Matrix<double, 6, 6> compute_material_tangent(const Matrix3d &F, int mat_index) override
+    Matrix<double, 6, 6> compute_material_tangent(const Matrix3d &F, int mat_index, ptrdiff_t element_idx, int i) override
     {
         Matrix3d C = compute_C(F);
         double   J = F.determinant();
@@ -91,16 +91,15 @@ class CompressibleNeoHookean : public LargeStrainMechModel {
         return C_tangent;
     }
 
-    Matrix<double, 9, 9> get_reference_stiffness() const override
+    Matrix<double, 9, 9> get_reference_stiffness() override
     {
         // Compute reference tangent at F=I
-        auto                *self       = const_cast<CompressibleNeoHookean *>(this);
         Matrix<double, 9, 9> kapparef   = Matrix<double, 9, 9>::Zero();
         Matrix3d             F_identity = Matrix3d::Identity();
 
         for (int mat_idx = 0; mat_idx < n_mat; ++mat_idx) {
-            Matrix3d             S        = self->compute_S(F_identity, mat_idx, 0);
-            Matrix<double, 6, 6> C_mandel = self->compute_material_tangent(F_identity, mat_idx);
+            Matrix3d             S        = compute_S(F_identity, mat_idx, 0, 0);
+            Matrix<double, 6, 6> C_mandel = compute_material_tangent(F_identity, mat_idx, 0, 0);
             Matrix<double, 9, 9> A        = compute_spatial_tangent(F_identity, S, C_mandel);
             kapparef += A;
         }