add J2 plasticity in hooke

linalg/dual_tensor.hpp

@@ -0,0 +1,49 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#pragma once
+
+#include "tensor.hpp"
+
+namespace mfem
+{
+namespace internal
+{
+
+template <typename value_type, typename gradient_type, int n>
+MFEM_HOST_DEVICE auto get_value(const
+                                tensor<dual<value_type, gradient_type>, n>& arg)
+{
+   tensor<double, n> output{};
+   for (int i = 0; i < n; i++)
+   {
+      output[i] = arg[i].value;
+   }
+   return output;
+}
+
+template <typename value_type, typename gradient_type, int m, int n>
+MFEM_HOST_DEVICE auto get_value(const
+                                tensor<dual<value_type, gradient_type>, m, n>& arg)
+{
+   tensor<double, m, n> output{};
+   for (int i = 0; i < m; i++)
+   {
+      for (int j = 0; j < n; j++)
+      {
+         output[i][j] = arg[i][j].value;
+      }
+   }
+   return output;
+}
+
+}
+}

linalg/tensor.hpp

@@ -1211,7 +1211,8 @@ T sqnorm(const tensor<T, m, n>& A)
 template <typename T, int... n> MFEM_HOST_DEVICE
 T norm(const tensor<T, n...>& A)
 {
-   return std::sqrt(sqnorm(A));
+   using std::sqrt;
+   return sqrt(sqnorm(A));
 }
 
 /**

miniapps/hooke/hooke.cpp

@@ -34,6 +34,7 @@
 
 #include "materials/linear_elastic.hpp"
 #include "materials/neohookean.hpp"
+#include "materials/j2_plastic.hpp"
 #include "operators/elasticity_gradient_operator.hpp"
 #include "operators/elasticity_operator.hpp"
 #include "preconditioners/diagonal_preconditioner.hpp"
@@ -128,8 +129,27 @@ int main(int argc, char *argv[])
    // * EnzymeRev
    // * FiniteDiff
    // * InternalFwd
-   const NeoHookeanMaterial<dimension, GradientType::InternalFwd> material{};
-   elasticity_op.SetMaterial(material);
+   // const NeoHookeanMaterial<dimension, GradientType::InternalFwd> material{};
+   double G = 79000;
+   double K = 10 * G;
+   double E  = 9 * K * G / (3 * K + G);
+   double nu = (3 * K - 2 * G) / (2 * (3 * K + G));
+   J2Material material
+   {
+      E,                // Young's modulus
+      nu,               // Poisson's ratio
+      10.0,              // isotropic hardening constant
+      0.0,              // kinematic hardening constant
+      165 * sqrt(3.0),  // yield stress
+      1.0               // mass density
+   };
+   auto material_state = elasticity_op.CreateMaterialState(material);
+   for (int i = 0; i < material_state.Size(); i++)
+   {
+      material_state[i] = {};
+   }
+
+   elasticity_op.SetMaterial(material, material_state);
 
    // Define all essential boundaries. In this specific example, this includes
    // all fixed and statically displaced degrees of freedom on mesh entities in
@@ -160,15 +180,12 @@ int main(int argc, char *argv[])
    Vector U;
    U_gf.GetTrueDofs(U);
 
-   // Prescribe a fixed displacement to the displaced degrees of freedom.
-   U.SetSubVector(elasticity_op.GetPrescribedDisplacementTDofs(), 1.0e-2);
-
    // Define the type of preconditioner to use for the linear solver.
    ElasticityDiagonalPreconditioner diagonal_pc(
       static_cast<ElasticityDiagonalPreconditioner::Type>(diagpc_type));
 
    CGSolver cg(MPI_COMM_WORLD);
-   cg.SetRelTol(1e-1);
+   cg.SetRelTol(1e-4);
    cg.SetMaxIter(10000);
    cg.SetPrintLevel(2);
    cg.SetPreconditioner(diagonal_pc);
@@ -176,14 +193,36 @@ int main(int argc, char *argv[])
    NewtonSolver newton(MPI_COMM_WORLD);
    newton.SetSolver(cg);
    newton.SetOperator(elasticity_op);
-   newton.SetRelTol(1e-6);
-   newton.SetMaxIter(10);
+   newton.SetRelTol(1e-8);
+   newton.SetAbsTol(1e-11);
+   newton.SetMaxIter(1000);
    newton.SetPrintLevel(1);
 
-   Vector zero;
-   newton.Mult(zero, U);
+   Array<double> displacement_increments(20);
 
-   U_gf.Distribute(U);
+   for (int i = 0; i < displacement_increments.Size(); i++)
+   {
+      printf("\n>>> solving displacement increment %d\n", i);
+      // Prescribe a fixed displacement to the displaced degrees of freedom.
+      U.SetSubVector(elasticity_op.GetPrescribedDisplacementTDofs(),
+                     min(i, displacement_increments.Size()-i-1)*1.0e-3);
+
+      Vector zero;
+      newton.Mult(zero, U);
+
+      elasticity_op.UpdateMaterialState();
+
+      for (int i = 0; i < material_state.Size(); i++)
+      {
+         if (material_state[i].accumulated_plastic_strain > 0.0)
+         {
+            printf(">>> plastic deformation reached\n");
+            break;
+         }
+      }
+
+      U_gf.Distribute(U);
+   }
 
    if (visualization)
    {

miniapps/hooke/kernels/elasticity_kernels.hpp

@@ -14,6 +14,7 @@
 
 #include "kernel_helpers.hpp"
 #include "linalg/vector.hpp"
+#include "../materials/material_traits.hpp"
 
 using mfem::internal::tensor;
 using mfem::internal::make_tensor;
@@ -23,7 +24,6 @@ namespace mfem
 
 namespace ElasticityKernels
 {
-
 /**
  * @brief Apply the 3D elasticity kernel
  *
@@ -45,7 +45,8 @@ namespace ElasticityKernels
  * @param Y_ Output vector d1d x d1d x d1d x vdim x ne.
  * @param material Material object.
  */
-template <int d1d, int q1d, typename material_type> static inline
+template <int d1d, int q1d, typename material_type, typename material_type_state>
+static inline
 void Apply3D(const int ne,
              const Array<double> &B_,
              const Array<double> &G_,
@@ -53,7 +54,8 @@ void Apply3D(const int ne,
              const Vector &Jacobian_,
              const Vector &detJ_,
              const Vector &X_, Vector &Y_,
-             const material_type &material)
+             const material_type &material,
+             Array<material_type_state> &material_state_)
 {
    static constexpr int dim = 3;
    KernelHelpers::CheckMemoryRestriction(d1d, q1d);
@@ -68,6 +70,8 @@ void Apply3D(const int ne,
    const auto J = Reshape(Jacobian_.Read(), q1d, q1d, q1d, dim, dim, ne);
    const auto detJ = Reshape(detJ_.Read(), q1d, q1d, q1d, ne);
    const auto U = Reshape(X_.Read(), d1d, d1d, d1d, dim, ne);
+   // This will access random memory if material_type_state is NoState
+   auto state = Reshape(material_state_.ReadWrite(), q1d, q1d, q1d, ne);
    auto force = Reshape(Y_.ReadWrite(), d1d, d1d, d1d, dim, ne);
 
    MFEM_FORALL_3D(e, ne, q1d, q1d, q1d,
@@ -95,10 +99,20 @@ void Apply3D(const int ne,
 
                // This represents the quadrature function operation in the
                // Finite Element Operator Decomposition e.g. A_Q(x).
-               auto sigma = material.stress(dudx);
+               if constexpr (material_has_state<material_type>::value)
+               {
+                  auto sigma = material.stress(state(qx,qy,qz,e), dudx);
 
-               invJ_sigma_detJw(qx, qy, qz) =
-                  invJqp * sigma * detJ(qx, qy, qz, e) * qweights(qx, qy, qz);
+                  invJ_sigma_detJw(qx, qy, qz) =
+                     invJqp * sigma * detJ(qx, qy, qz, e) * qweights(qx, qy, qz);
+               }
+               else
+               {
+                  auto sigma = material.stress(dudx);
+
+                  invJ_sigma_detJw(qx, qy, qz) =
+                     invJqp * sigma * detJ(qx, qy, qz, e) * qweights(qx, qy, qz);
+               }
             }
          }
       }
@@ -138,12 +152,14 @@ void Apply3D(const int ne,
  * @param dsigma_cache_ Vector to use as memory for the cache of dsigma/dudx.
  * Size needed is ne x q1d x q1d x q1d x dim x dim x dim x dim.
  */
-template <int d1d, int q1d, typename material_type> static inline
+template <int d1d, int q1d, typename material_type, typename material_type_state>
+static inline
 void ApplyGradient3D(const int ne,
                      const Array<double> &B_, const Array<double> &G_,
                      const Array<double> &W_, const Vector &Jacobian_,
                      const Vector &detJ_, const Vector &dU_, Vector &dF_,
                      const Vector &U_, const material_type &material,
+                     Array<material_type_state> &material_state_,
                      const bool use_cache_, const bool recompute_cache_,
                      Vector &dsigma_cache_)
 {
@@ -162,6 +178,8 @@ void ApplyGradient3D(const int ne,
    const auto dU = Reshape(dU_.Read(), d1d, d1d, d1d, dim, ne);
    auto force = Reshape(dF_.ReadWrite(), d1d, d1d, d1d, dim, ne);
    const auto U = Reshape(U_.Read(), d1d, d1d, d1d, dim, ne);
+   // This will access random memory if material_type_state is NoState
+   auto state = Reshape(material_state_.ReadWrite(), q1d, q1d, q1d, ne);
 
    auto dsigma_cache = Reshape(dsigma_cache_.ReadWrite(), ne, q1d, q1d, q1d,
                                dim, dim, dim, dim);
@@ -204,7 +222,7 @@ void ApplyGradient3D(const int ne,
 
                   if (recompute_cache_)
                   {
-                     C = material.gradient(dudx);
+                     C = material.gradient(state(qx,qy,qz,e), dudx);
                      for (int i = 0; i < dim; i++)
                      {
                         for (int j = 0; j < dim; j++)
@@ -229,7 +247,8 @@ void ApplyGradient3D(const int ne,
                }
                else
                {
-                  auto dsigma = material.action_of_gradient(dudx, ddudx);
+                  auto dsigma = material.action_of_gradient(state(qx,qy,qz,e), dudx,
+                                                            ddudx);
                   invJ_dsigma_detJw(qx, qy, qz) =
                      invJqp * dsigma * detJ(qx, qy, qz, e) * qweights(qx, qy, qz);
                }
@@ -265,7 +284,8 @@ void ApplyGradient3D(const int ne,
  * matrices. Size needed is d1d x d1d x d1d x dim x ne x dim.
  * @param material Material object.
  */
-template <int d1d, int q1d, typename material_type> static inline
+template <int d1d, int q1d, typename material_type, typename material_type_state>
+static inline
 void AssembleGradientDiagonal3D(const int ne,
                                 const Array<double> &B_,
                                 const Array<double> &G_,
@@ -274,7 +294,8 @@ void AssembleGradientDiagonal3D(const int ne,
                                 const Vector &detJ_,
                                 const Vector &X_,
                                 Vector &Ke_diag_memory,
-                                const material_type &material)
+                                const material_type &material,
+                                Array<material_type_state> &material_state)
 {
    static constexpr int dim = 3;
    KernelHelpers::CheckMemoryRestriction(d1d, q1d);
@@ -289,6 +310,8 @@ void AssembleGradientDiagonal3D(const int ne,
    const auto J = Reshape(Jacobian_.Read(), q1d, q1d, q1d, dim, dim, ne);
    const auto detJ = Reshape(detJ_.Read(), q1d, q1d, q1d, ne);
    const auto U = Reshape(X_.Read(), d1d, d1d, d1d, dim, ne);
+   // This will access random memory if material_type_state is NoState
+   auto state = Reshape(material_state.ReadWrite(), q1d, q1d, q1d, ne);
 
    auto Ke_diag_m =
       Reshape(Ke_diag_memory.ReadWrite(), d1d, d1d, d1d, dim, ne, dim);
@@ -318,10 +341,11 @@ void AssembleGradientDiagonal3D(const int ne,
 
                const auto dudx = dudxi(qz, qy, qx) * invJqp;
 
-               const auto dsigma_ddudx = material.gradient(dudx);
+               const auto dsigma_ddudx = material.gradient(state(qx,qy,qz,e), dudx);
 
                const double JxW = detJ(qx, qy, qz, e) * qweights(qx, qy, qz);
-               const auto dphidx = KernelHelpers::GradAllShapeFunctions(qx, qy, qz, B, G, invJqp);
+               const auto dphidx = KernelHelpers::GradAllShapeFunctions(qx, qy, qz, B, G,
+                                                                        invJqp);
 
                for (int dx = 0; dx < d1d; dx++)
                {

miniapps/hooke/materials/j2_plastic.hpp

@@ -0,0 +1,181 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#pragma once
+
+#include "linalg/dual_tensor.hpp"
+#include "material_traits.hpp"
+#include "mfem.hpp"
+
+struct J2Material
+{
+   /// this material is written for 3D
+   static constexpr int dim = 3;
+
+   double E;        ///< Young's modulus
+   double nu;       ///< Poisson's ratio
+   double Hi;       ///< isotropic hardening constant
+   double Hk;       ///< kinematic hardening constant
+   double sigma_y;  ///< yield stress
+   double density;  ///< mass density
+
+   /// @brief variables required to characterize the hysteresis response
+   struct State
+   {
+      /// back-stress tensor
+      mfem::internal::tensor<double, dim, dim> beta;
+
+      /// plastic strain
+      mfem::internal::tensor<double, dim, dim> plastic_strain;
+
+      /// accumulated plastic strain
+      double accumulated_plastic_strain;
+   };
+
+   /** @brief calculate the Cauchy stress, given the displacement gradient and previous material state */
+   template <typename T>
+   auto stress(State& state, const T du_dx) const
+   {
+      using std::sqrt;
+      constexpr auto I = mfem::internal::IsotropicIdentity<3>();
+      const double   K = E / (3.0 * (1.0 - 2.0 * nu));
+      const double   G = 0.5 * E / (1.0 + nu);
+
+      //
+      // see pg. 260, box 7.5,
+      // in "Computational Methods for Plasticity"
+      //
+
+      // (i) elastic predictor
+      auto el_strain = mfem::internal::sym(du_dx) - state.plastic_strain;
+      auto p         = K * mfem::internal::tr(el_strain);
+      auto s         = 2.0 * G * mfem::internal::dev(el_strain);
+      auto eta       = s - state.beta;
+      auto q         = sqrt(3.0 / 2.0) * mfem::internal::norm(eta);
+      auto phi       = q - (sigma_y + Hi * state.accumulated_plastic_strain);
+
+      // (ii) admissibility
+      if (phi > 0.0)
+      {
+         // see (7.207) on pg. 261
+         auto plastic_strain_inc = phi / (3 * G + Hk + Hi);
+
+         // from here on, only normalize(eta) is required
+         // so we overwrite eta with its normalized version
+         eta = mfem::internal::normalize(eta);
+
+         // (iii) return mapping
+         s = s - sqrt(6.0) * G * plastic_strain_inc * eta;
+         state.accumulated_plastic_strain += mfem::internal::get_value(
+                                                plastic_strain_inc);
+         state.plastic_strain += sqrt(3.0 / 2.0) * mfem::internal::get_value(
+                                    plastic_strain_inc) *
+                                 mfem::internal::get_value(eta);
+         state.beta = state.beta + sqrt(2.0 / 3.0) * Hk * mfem::internal::get_value(
+                         plastic_strain_inc) * mfem::internal::get_value(eta);
+      }
+
+      return s + p * I;
+   }
+
+   /**
+    * @brief Apply the gradient of the stress.
+    *
+    */
+   MFEM_HOST_DEVICE
+   template <typename T> auto
+   action_of_gradient(State &state, const T &dudx, const T &ddudx) const
+   {
+      auto sigma = stress(state, make_tensor<dim, dim>([&](int i, int j)
+      {
+         return mfem::internal::dual<double, double> {dudx[i][j], ddudx[i][j]};
+      }));
+      return make_tensor<dim, dim>(
+      [&](int i, int j) { return sigma[i][j].gradient; });
+   }
+
+   MFEM_HOST_DEVICE tensor<double, dim, dim, dim, dim>
+   gradient(State &state, tensor<double, dim, dim> dudx) const
+   {
+      using std::sqrt;
+      constexpr auto I = mfem::internal::IsotropicIdentity<3>();
+      const double   K = E / (3.0 * (1.0 - 2.0 * nu));
+      const double   G = 0.5 * E / (1.0 + nu);
+
+      //
+      // see pg. 260, box 7.5,
+      // in "Computational Methods for Plasticity"
+      //
+
+      // (i) elastic predictor
+      auto el_strain = mfem::internal::sym(dudx) - state.plastic_strain;
+      auto p         = K * mfem::internal::tr(el_strain);
+      auto s         = 2.0 * G * mfem::internal::dev(el_strain);
+      auto eta       = s - state.beta;
+      auto q         = sqrt(3.0 / 2.0) * mfem::internal::norm(eta);
+      auto phi       = q - (sigma_y + Hi * state.accumulated_plastic_strain);
+
+      // (ii) admissibility
+      if (phi > 0.0)
+      {
+         // see (7.207) on pg. 261
+         auto plastic_strain_inc = phi / (3 * G + Hk + Hi);
+
+         // from here on, only normalize(eta) is required
+         // so we overwrite eta with its normalized version
+         eta = mfem::internal::normalize(eta);
+
+         // (iii) return mapping
+         s = s - sqrt(6.0) * G * plastic_strain_inc * eta;
+         state.accumulated_plastic_strain += mfem::internal::get_value(
+                                                plastic_strain_inc);
+         state.plastic_strain += sqrt(3.0 / 2.0) * mfem::internal::get_value(
+                                    plastic_strain_inc) *
+                                 mfem::internal::get_value(eta);
+         state.beta = state.beta + sqrt(2.0 / 3.0) * Hk * mfem::internal::get_value(
+                         plastic_strain_inc) * mfem::internal::get_value(eta);
+      }
+
+      // return s + p * I;
+
+      double A1 = 2.0 * G;
+      double A2 = 0.0;
+
+      tensor<double, 3, 3> N{};
+
+      auto pl_strain_inc = fmax(0.0, phi / (3 * G + Hk + Hi));
+
+      if (pl_strain_inc > 0.0)
+      {
+         tensor<double, 3, 3> s = 2.0 * G * mfem::internal::dev(el_strain);
+         N = mfem::internal::normalize(s - state.beta);
+
+         A1 -= 6 * G * G * pl_strain_inc / q;
+         A2 = 6 * G * G * ((pl_strain_inc / q) - (1.0 /
+                                                  (3.0 * G + Hi + Hk)));
+      }
+
+      return make_tensor<3, 3, 3, 3>([&](auto i, auto j, auto k, auto l)
+      {
+         double I4    = (i == j) * (k == l);
+         double I4sym = 0.5 * ((i == k) * (j == l) + (i == l) * (j == k));
+         double I4dev = I4sym - (i == j) * (k == l) / 3.0;
+         return K * I4 + A1 * I4dev + A2 * N(i, j) * N(k, l);
+      });
+   }
+
+};
+
+template<>
+struct material_has_state<J2Material>
+{
+   static const bool value = true;
+};

miniapps/hooke/materials/material_traits.hpp

@@ -0,0 +1,24 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#pragma once
+
+#include <mfem.hpp>
+
+struct NoState {};
+// <3 cpp...
+inline mfem::Array<NoState> NoStateArray;
+
+template<typename T>
+struct material_has_state
+{
+   static const bool value = false;
+};

miniapps/hooke/operators/elasticity_operator.cpp

@@ -86,7 +86,7 @@ void ElasticityOperator::Mult(const Vector &X, Vector &Y) const
    // Apply operator
    element_apply_kernel_wrapper(ne_, maps_->B, maps_->G, ir_->GetWeights(),
                                 geometric_factors_->J, geometric_factors_->detJ,
-                                X_el_, Y_el_);
+                                X_el_, Y_el_, false);
 
    // E-vector to L-vector
    h1_element_restriction_->MultTranspose(Y_el_, Y_local_);

miniapps/hooke/operators/elasticity_operator.hpp

@@ -13,6 +13,7 @@
 #define MFEM_ELASTICITY_OP_HPP
 
 #include "../kernels/elasticity_kernels.hpp"
+#include "../materials/material_traits.hpp"
 #include "mfem.hpp"
 
 namespace mfem
@@ -83,6 +84,12 @@ public:
    void AssembleGradientDiagonal(Vector &Ke_diag, Vector &K_diag_local,
                                  Vector &K_diag) const;
 
+   template<typename material_type>
+   Array<typename material_type::State> CreateMaterialState(material_type material)
+   {
+      return Array<typename material_type::State>(q1d_*q1d_*q1d_*ne_);
+   }
+
    ~ElasticityOperator();
 
    ParMesh &mesh_;
@@ -143,7 +150,7 @@ public:
     */
    std::function<void(const int, const Array<double> &, const Array<double> &,
                       const Array<double> &, const Vector &, const Vector &,
-                      const Vector &, Vector &)>
+                      const Vector &, Vector &, bool retrieve_swap)>
    element_apply_kernel_wrapper;
 
    /**
@@ -175,8 +182,9 @@ public:
     * @tparam material_type
     * @param[in] material
     */
-   template <typename material_type>
-   void SetMaterial(const material_type &material)
+   template <typename material_type, typename material_type_state = NoState>
+   void SetMaterial(const material_type &material,
+                    Array<material_type_state> &material_state = NoStateArray)
    {
       if (dim_ != 3)
       {
@@ -184,89 +192,120 @@ public:
       }
 
       element_apply_kernel_wrapper =
-         [=](const int ne, const Array<double> &B_, const Array<double> &G_,
-             const Array<double> &W_, const Vector &Jacobian_,
-             const Vector &detJ_, const Vector &X_, Vector &Y_)
+         [=, material_state_swap = material_state, &material_state]
+         (const int ne,
+          const Array<double> &B_,
+          const Array<double> &G_,
+          const Array<double> &W_, const Vector &Jacobian_,
+          const Vector &detJ_, const Vector &X_, Vector &Y_, bool retrieve_swap = false)
+         mutable
       {
+         if (retrieve_swap)
+         {
+            material_state = material_state_swap;
+            return;
+         }
+         material_state_swap = material_state;
+
          const int id = (d1d_ << 4) | q1d_;
          switch (id)
          {
             case 0x22:
             {
                ElasticityKernels::Apply3D<2, 2, material_type>(
-                  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
+                  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material, material_state_swap);
                break;
             }
-            case 0x33:
-            {
-               ElasticityKernels::Apply3D<3, 3, material_type>(
-                  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
-               break;
-            }
-            case 0x44:
-               ElasticityKernels::Apply3D<4, 4, material_type>(
-                  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
-               break;
+            // case 0x33:
+            // {
+            //    ElasticityKernels::Apply3D<3, 3, material_type>(
+            //       ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material, material_state_swap);
+            //    break;
+            // }
+            // case 0x44:
+            //    ElasticityKernels::Apply3D<4, 4, material_type>(
+            //       ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material, material_state_swap);
+            //    break;
             default:
                MFEM_ABORT("Not implemented: " << std::hex << id << std::dec);
          }
       };
 
       element_apply_gradient_kernel_wrapper =
-         [=](const int ne, const Array<double> &B_, const Array<double> &G_,
-             const Array<double> &W_, const Vector &Jacobian_,
-             const Vector &detJ_, const Vector &dU_, Vector &dF_,
-             const Vector &U_)
+         [=, material_state_swap = material_state, &material_state]
+         (const int ne,
+          const Array<double> &B_,
+          const Array<double> &G_,
+          const Array<double> &W_, const Vector &Jacobian_,
+          const Vector &detJ_, const Vector &dU_, Vector &dF_,
+          const Vector &U_)
+         mutable
       {
+         material_state_swap = material_state;
+
          const int id = (d1d_ << 4) | q1d_;
          switch (id)
          {
             case 0x22:
                ElasticityKernels::ApplyGradient3D<2, 2, material_type>(
-                  ne, B_, G_, W_, Jacobian_, detJ_, dU_, dF_, U_, material,
-                  use_cache_, recompute_cache_, dsigma_cache_);
-               break;
-            case 0x33:
-               ElasticityKernels::ApplyGradient3D<3, 3, material_type>(
-                  ne, B_, G_, W_, Jacobian_, detJ_, dU_, dF_, U_, material,
-                  use_cache_, recompute_cache_, dsigma_cache_);
-               break;
-            case 0x44:
-               ElasticityKernels::ApplyGradient3D<4, 4, material_type>(
-                  ne, B_, G_, W_, Jacobian_, detJ_, dU_, dF_, U_, material,
+                  ne, B_, G_, W_, Jacobian_, detJ_, dU_, dF_, U_, material, material_state_swap,
                   use_cache_, recompute_cache_, dsigma_cache_);
                break;
+            // case 0x33:
+            //    ElasticityKernels::ApplyGradient3D<3, 3, material_type>(
+            //       ne, B_, G_, W_, Jacobian_, detJ_, dU_, dF_, U_, material,
+            //       use_cache_, recompute_cache_, dsigma_cache_);
+            //    break;
+            // case 0x44:
+            //    ElasticityKernels::ApplyGradient3D<4, 4, material_type>(
+            //       ne, B_, G_, W_, Jacobian_, detJ_, dU_, dF_, U_, material,
+            //       use_cache_, recompute_cache_, dsigma_cache_);
+            //    break;
             default:
                MFEM_ABORT("Not implemented for D1D=" << d1d_ << " and Q1D=" << q1d_);
          }
       };
 
       element_kernel_assemble_diagonal_wrapper =
-         [=](const int ne, const Array<double> &B_, const Array<double> &G_,
-             const Array<double> &W_, const Vector &Jacobian_,
-             const Vector &detJ_, const Vector &X_, Vector &Y_)
+         [=, material_state_swap = material_state, &material_state]
+         (const int ne, const Array<double> &B_,
+          const Array<double> &G_,
+          const Array<double> &W_, const Vector &Jacobian_,
+          const Vector &detJ_, const Vector &X_, Vector &Y_)
+         mutable
       {
+         material_state_swap = material_state;
+
          const int id = (d1d_ << 4) | q1d_;
          switch (id)
          {
             case 0x22:
                ElasticityKernels::AssembleGradientDiagonal3D<2, 2, material_type>(
-                  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
-               break;
-            case 0x33:
-               ElasticityKernels::AssembleGradientDiagonal3D<3, 3, material_type>(
-                  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
-               break;
-            case 0x44:
-               ElasticityKernels::AssembleGradientDiagonal3D<4, 4, material_type>(
-                  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
+                  ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material, material_state_swap);
                break;
+            // case 0x33:
+            //    ElasticityKernels::AssembleGradientDiagonal3D<3, 3, material_type>(
+            //       ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
+            //    break;
+            // case 0x44:
+            //    ElasticityKernels::AssembleGradientDiagonal3D<4, 4, material_type>(
+            //       ne, B_, G_, W_, Jacobian_, detJ_, X_, Y_, material);
+            //    break;
             default:
                MFEM_ABORT("Not implemented: " << std::hex << id << std::dec);
          }
       };
    }
 
+   void UpdateMaterialState()
+   {
+      // This is a dummy call. We only retrieve the latest material state and
+      // exit.
+      element_apply_kernel_wrapper(ne_, maps_->B, maps_->G, ir_->GetWeights(),
+                                   geometric_factors_->J,
+                                   geometric_factors_->detJ, X_el_, Y_el_, true);
+   }
+
    /**
     * @brief Set the essential attributes which mark degrees of freedom for the
     * solving process.