diff --git a/src/model/materials/internal.hh b/src/model/materials/internal.hh
index e71461b..8b93a7d 100644
--- a/src/model/materials/internal.hh
+++ b/src/model/materials/internal.hh
@@ -1,63 +1,63 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef INTERNAL_HH
 #define INTERNAL_HH
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include <memory>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 template <typename T, UInt dim>
 struct Internal {
   /// Stored field
   std::shared_ptr<Grid<T, dim>> field;
 
   /// Initialize previous field
   void initialize() {
     previous = std::make_shared<Grid<T, dim>>(field->sizes(),
                                               field->getNbComponents());
     update();
   }
 
   /// Update the field previous state
   void update() { *previous = *field; }
   /// Reset the current field
   void reset() { *field = *previous; }
 
   /// Assign the field pointer
   void operator=(decltype(field)&& field) { this->field = std::move(field); }
   /// Dereference the field pointer
   typename decltype(field)::element_type& operator*() { return *field; }
   /// Dereference the field pointer (const version)
   const typename decltype(field)::element_type& operator*() const {
     return *field;
   }
   /// Upcast to GridBase pointer
   operator std::shared_ptr<GridBase<T>>() { return field; }
 
-private:
+public:
   decltype(field) previous;
 };
 
 }  // namespace tamaas
 #endif  // INTERNAL_HH
diff --git a/src/model/materials/isotropic_hardening.cpp b/src/model/materials/isotropic_hardening.cpp
index e6df614..62a5e66 100644
--- a/src/model/materials/isotropic_hardening.cpp
+++ b/src/model/materials/isotropic_hardening.cpp
@@ -1,178 +1,195 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "isotropic_hardening.hh"
 #include "influence.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 IsotropicHardening<type>::IsotropicHardening(Model* model, Real sigma_0, Real h)
     : Material<type>(model), sigma_0(sigma_0), h(h) {
   plastic_strain = model->request<type, false, Real>(
       "plastic_strain", model->getDiscretization(), trait::voigt);
   cumulated_plastic_strain = model->request<type, false, Real>(
       "cumulated_plastic_strain", model->getDiscretization(), 1);
 
   plastic_strain.initialize();
   cumulated_plastic_strain.initialize();
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void IsotropicHardening<model_type::volume_2d>::
     computeInelasticDeformationIncrement(Field& increment, const Field& strain,
                                          const Field& strain_increment) {
   influence::ElasticHelper<dim> elasticity(this->model->getShearModulus(),
                                            this->model->getPoissonRatio());
   plastic_strain.reset();
   cumulated_plastic_strain.reset();
 
   Real sigma_0 = this->sigma_0, h = this->h;  // for captures
 
   Loop::loop(
       [elasticity, sigma_0,
        h] CUDA_LAMBDA(Mat dep,             // < plastic strain increment
                       CMat epsilon,        // < total strain
                       CMat delta_epsilon,  // < strain increment
                       Mat ep,              // < plastic strain
                       Real & p /* < cumulated plastic strain */) {
         // Trial elastic stress
         auto sigma_tr = elasticity(epsilon - ep + delta_epsilon);
 
         // Zero increment by default
         dep = 0;
 
         // Deviatoric trial stress
         decltype(sigma_tr) dev;
         dev.deviatoric(sigma_tr);
 
         // Von-Mises stress
         Real von_mises = std::sqrt(1.5) * dev.l2norm();
 
         // Plasticity condition
         if (von_mises - IsotropicHardening<model_type::volume_2d>::hardening(
                             p, h, sigma_0) >
             0) {
           // radial return
           Real dp =
               (von_mises - IsotropicHardening<model_type::volume_2d>::hardening(
                                p, h, sigma_0)) /
               (3 * elasticity.mu + h);
 
           dev *= 3 * dp / (2 * von_mises);  // saving memory (dev is delta ep)
           dep = dev;
 
           p += dp;
           ep += dep;
         }
       },
       range<Mat>(increment), range<CMat>(strain), range<CMat>(strain_increment),
       range<Mat>(*plastic_strain), *cumulated_plastic_strain);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void IsotropicHardening<model_type::volume_2d>::computeStress(
     Field& stress, const Field& strain, const Field& strain_increment) {
   this->computeInelasticDeformationIncrement(stress, strain, strain_increment);
 
   const influence::ElasticHelper<dim> elasticity(
       this->model->getShearModulus(), this->model->getPoissonRatio());
 
   Loop::loop(
       [elasticity] CUDA_LAMBDA(Mat sigma,           // < stress
                                CMat epsilon,        // < total strain
                                CMat delta_epsilon,  // < strain increment
                                CMat ep              // < plastic strain
       ) {
         sigma = epsilon;
         sigma += delta_epsilon;
         sigma -= ep;
         sigma = elasticity(sigma);
       },
       range<Mat>(stress), range<CMat>(strain), range<CMat>(strain_increment),
       range<CMat>(*plastic_strain));
 }
 
+/* -------------------------------------------------------------------------- */
+template <>
+void IsotropicHardening<model_type::volume_2d>::computeEigenStress(
+    Field& stress, const Field& strain, const Field& strain_increment) {
+  this->computeInelasticDeformationIncrement(stress, strain, strain_increment);
+
+  const influence::ElasticHelper<dim> elasticity(
+      this->model->getShearModulus(), this->model->getPoissonRatio());
+
+  Loop::loop([elasticity] CUDA_LAMBDA(Mat sigma,    // < stress
+                                      CMat ep,      // < plastic strain
+                                      CMat ep_prev  // < previous plastic strain
+             ) { sigma = elasticity(ep - ep_prev); },
+             range<Mat>(stress), range<CMat>(*plastic_strain),
+             range<CMat>(*plastic_strain.previous));
+}
+
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void IsotropicHardening<type>::update() {
   plastic_strain.update();
   cumulated_plastic_strain.update();
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void IsotropicHardening<model_type::volume_2d>::applyTangentIncrement(
     Grid<Real, dim>& output, const Grid<Real, dim>& input,
     const Grid<Real, dim>& strain,
     const Grid<Real, dim>& strain_increment) const {
   const influence::ElasticHelper<dim> elasticity(
       this->model->getShearModulus(), this->model->getPoissonRatio());
 
   Real sigma_0 = this->sigma_0, h = this->h;  // for captures
 
   Loop::loop(
       [elasticity, sigma_0, h] CUDA_LAMBDA(Mat out, CMat in, CMat epsilon,
                                            CMat delta_epsilon, CMat ep,
                                            const Real& p) {
         auto sigma_tr = elasticity(epsilon - ep + delta_epsilon);
 
         decltype(sigma_tr) dev;
         dev.deviatoric(sigma_tr, 3);
 
         auto von_mises = std::sqrt(1.5) * dev.l2norm();
 
         if (von_mises - IsotropicHardening<model_type::volume_2d>::hardening(
                             p, h, sigma_0) >
             0) {
           Real dp =
               (von_mises - IsotropicHardening<model_type::volume_2d>::hardening(
                                p, h, sigma_0)) /
               (3 * elasticity.mu + h);
           // Applying tangent from Bonnet & Frangi, p.175
 
           const Real beta = 3 * elasticity.mu * dp / von_mises;
           const Real gamma = 3 * elasticity.mu / (3 * elasticity.mu + h);
 
           const Real dot = dev.dot(in);
           dev *= 3 * elasticity.mu * (gamma - beta) * dot /
                  (von_mises * von_mises);
 
           out.deviatoric(in);
           out *= 2 * elasticity.mu * beta;
           out += dev;
         } else
           out = 0;
       },
       range<Mat>(output), range<CMat>(input), range<CMat>(strain),
       range<CMat>(strain_increment), range<CMat>(*this->plastic_strain),
       *this->cumulated_plastic_strain);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation */
 /* -------------------------------------------------------------------------- */
 template class IsotropicHardening<model_type::volume_2d>;
 
 }  // namespace tamaas
diff --git a/src/model/materials/isotropic_hardening.hh b/src/model/materials/isotropic_hardening.hh
index 3191134..5d0253e 100644
--- a/src/model/materials/isotropic_hardening.hh
+++ b/src/model/materials/isotropic_hardening.hh
@@ -1,90 +1,94 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef ISOTROPIC_HARDENING_HH
 #define ISOTROPIC_HARDENING_HH
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "influence.hh"
 #include "internal.hh"
 #include "material.hh"
 #include "model.hh"
 #include "model_type.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 class IsotropicHardening : Material<type> {
   using parent = Material<type>;
   using trait = typename parent::trait;
   static constexpr UInt dim = parent::trait::dimension;
   using Field = typename parent::Field;
 
 public:
   /// Constructor
   IsotropicHardening(Model* model, Real sigma_0, Real h);
 
   /// Compute plastic strain increment with radial return algorithm
   void
   computeInelasticDeformationIncrement(Field& increment, const Field& strain,
                                        const Field& strain_increment) override;
 
   /// Compute stress
   void computeStress(Field& stress, const Field& strain,
                      const Field& strain_increment) override;
 
+  /// Compute stress due to plastic strain increment
+  void computeEigenStress(Field& stress, const Field& strain,
+                          const Field& strain_increment) override;
+
   /// Update internal variables
   void update() override;
 
   void applyTangentIncrement(Field& output, const Field& input,
                              const Field& strain,
                              const Field& strain_increment) const;
 
   /// Linear hardening function
   static __device__ __host__ Real hardening(Real p, Real h, Real sigma_0) {
     return sigma_0 + h * p;
   }
 
   Real getHardeningModulus() const { return h; }
   Real getYieldStress() const { return sigma_0; }
   const GridBase<Real>& getPlasticStrain() const { return *plastic_strain; }
   GridBase<Real>& getPlasticStrain() { return *plastic_strain; }
 
   void setHardeningModulus(Real h_) {
     if (h_ < 0)
       TAMAAS_EXCEPTION("Hardening modulus should be positive");
     h = h_;
   }
   void setYieldStress(Real sigma_0_) {
     if (sigma_0_ < 0)
       TAMAAS_EXCEPTION("Yield stress should be positive");
     sigma_0 = sigma_0_;
   }
 
 protected:
   Real sigma_0;  /// < initial yield stress
   Real h;        /// < hardening modulus
   Internal<Real, 3> plastic_strain, cumulated_plastic_strain;
 };
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
 #endif  // ISOTROPIC_HARDENING_HH
diff --git a/src/model/materials/material.hh b/src/model/materials/material.hh
index af73116..67460f4 100644
--- a/src/model/materials/material.hh
+++ b/src/model/materials/material.hh
@@ -1,67 +1,71 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef MATERIAL_HH
 #define MATERIAL_HH
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "model.hh"
 #include "model_type.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 template <model_type type>
 class Material {
   // Defining some helper types
 protected:
   using trait = model_type_traits<type>;
   static constexpr UInt dim = trait::dimension;
   using Field = Grid<Real, dim>;
   using Mat = SymMatrixProxy<Real, dim>;
   using CMat = SymMatrixProxy<const Real, dim>;
 
 public:
   /// Constructor
   Material(Model* model) : model(model) {}
   /// Destructor
   virtual ~Material() = default;
 
   // Material Interface
 public:
   /// Compute the increment of inelastic deformation from the total strain and
   /// total strain increment
   virtual void
   computeInelasticDeformationIncrement(Field& increment, const Field& strain,
                                        const Field& strain_increment) = 0;
 
   /// Compute the stress from total strain and strain increment
   virtual void computeStress(Field& stress, const Field& strain,
                              const Field& strain_increment) = 0;
 
+  /// Compute stress due to inelastic increment
+  virtual void computeEigenStress(Field& stress, const Field& strain,
+                                  const Field& strain_increment) = 0;
+
   /// Update internal variables
   virtual void update() = 0;
 
 protected:
   Model* model;
 };
 
 }  // namespace tamaas
 #endif  // MATERIAL_HH
diff --git a/src/model/residual.cpp b/src/model/residual.cpp
index b7fcd4f..72a4f46 100644
--- a/src/model/residual.cpp
+++ b/src/model/residual.cpp
@@ -1,159 +1,158 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "residual.hh"
 #include "grid_view.hh"
 #include "model_factory.hh"
 #include "model_type.hh"
 #include <list>
 #include <memory>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 ResidualTemplate<type>::ResidualTemplate(Model* model, Real sigma_0, Real h)
     : Residual(model), hardening(model, sigma_0, h) {
 
   // Registering operators for residual and displacement calculation
   ModelFactory::registerVolumeOperators(*model);
 
   const auto n = model->getDiscretization();
   const auto nc = trait::voigt;
 
   strain = model->request<type, false, Real>("strain", n, nc);
   stress = model->request<type, false, Real>("stress", n, nc);
   residual = model->request<type, false, Real>("residual", n, nc);
   tmp = allocateGrid<type, false, Real>(n, nc);
 
   plastic_filter = [this](UInt layer) {
 #if 0
     return true;
 #else
     return this->plastic_layers.find(layer) != this->plastic_layers.end();
 #endif
   };
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::computeResidual(GridBase<Real>& strain_increment) {
   Grid<Real, dim> inc(strain->sizes(), strain->getNbComponents(),
                       strain_increment.view());
-  hardening.computeInelasticDeformationIncrement(*residual, *strain, inc);
+  hardening.computeEigenStress(*residual, *strain, inc);
   this->updateFilter(*residual);
-  this->model->applyElasticity(*residual, *residual);
   integralOperator("mindlin_gradient")
       .applyIf(*residual, *residual, plastic_filter);
   integralOperator("boussinesq_gradient")
       .apply(this->model->getTraction(), *tmp);
   *residual -= strain_increment;
   *residual += *tmp;
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::applyTangent(
     GridBase<Real>& output, GridBase<Real>& input,
     GridBase<Real>& current_strain_increment) {
   auto& inc = dynamic_cast<Grid<Real, dim>&>(current_strain_increment);
   auto& out = dynamic_cast<Grid<Real, dim>&>(output);
   auto& in = dynamic_cast<Grid<Real, dim>&>(input);
 
   hardening.applyTangentIncrement(out, in, *strain, inc);
   this->updateFilter(out);
   integralOperator("mindlin_gradient").applyIf(out, out, plastic_filter);
   out -= in;
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::computeStress(GridBase<Real>& strain_increment) {
   auto& inc = dynamic_cast<Grid<Real, dim>&>(strain_increment);
   hardening.computeStress(*stress, *strain, inc);
   // integralOperator("boussinesq_gradient")
   //     .apply(this->model->getTraction(), tmp);
   // this->model->applyElasticity(tmp, tmp);
   // stress += tmp;
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::updateState(
     GridBase<Real>& converged_strain_increment) {
   auto& inc = dynamic_cast<Grid<Real, dim>&>(converged_strain_increment);
   hardening.computeStress(*stress, *strain, inc);
   hardening.update();
   this->model->applyElasticity(*residual, hardening.getPlasticStrain());
   updateFilter(*residual);
   *strain += converged_strain_increment;
 
   // Computing displacements
   integralOperator("mindlin").applyIf(*residual, this->model->getDisplacement(),
                                       plastic_filter);
   Grid<Real, dim> disp_tmp(this->model->getDiscretization(), trait::components);
   integralOperator("boussinesq").apply(this->model->getTraction(), disp_tmp);
   this->model->getDisplacement() += disp_tmp;
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::computeResidualDisplacement(
     GridBase<Real>& strain_increment) {
   auto& inc = dynamic_cast<Grid<Real, dim>&>(strain_increment);
   hardening.computeInelasticDeformationIncrement(*tmp, *strain, inc);
   updateFilter(*tmp);
   this->model->applyElasticity(*residual, *tmp);
 
   integralOperator("mindlin").applyIf(*residual, this->model->getDisplacement(),
                                       plastic_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::updateFilter(
     Grid<Real, dim>& plastic_strain_increment) {
   // plastic_layers.clear();
   for (UInt i : Loop::range(plastic_strain_increment.sizes().front())) {
     auto slice = make_view(plastic_strain_increment, i);
     if (slice.dot(slice) / slice.globalDataSize() > 1e-14)
       plastic_layers.insert(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::setIntegrationMethod(integration_method method,
                                                   Real cutoff) {
   auto& mindlin = dynamic_cast<Mindlin<type, 1>&>(
       *this->model->getIntegralOperator("mindlin"));
   auto& mindlin_g = dynamic_cast<Mindlin<type, 2>&>(
       *this->model->getIntegralOperator("mindlin_gradient"));
 
   mindlin.setIntegrationMethod(method, cutoff);
   mindlin_g.setIntegrationMethod(method, cutoff);
 }
 
 /* -------------------------------------------------------------------------- */
 template class ResidualTemplate<model_type::volume_2d>;
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas