diff --git a/python/wrap/test_features.cpp b/python/wrap/test_features.cpp
index a8fba4b..6d4c0ae 100644
--- a/python/wrap/test_features.cpp
+++ b/python/wrap/test_features.cpp
@@ -1,67 +1,93 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "kelvin.hh"
 #include "mindlin.hh"
 #include "model.hh"
 #include "model_type.hh"
 #include "volume_potential.hh"
+#include "isotropic_hardening.hh"
 #include "wrap.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 namespace wrap {
 
 template <template <model_type, UInt> class KOp, model_type type,
           UInt tensor_order>
 void wrapKOp(const std::string& basename, py::module& mod) {
   constexpr UInt dim = model_type_traits<type>::dimension;
   std::stringstream str;
   str << basename << tensor_order;
   py::class_<KOp<type, tensor_order>>(mod, str.str().c_str())
       .def(py::init<Model*>())
       .def("apply",
            [](const KOp<type, tensor_order>& engine, Grid<Real, dim>& in,
               Grid<Real, dim>& out) { engine.apply(in, out); });
 }
 
+template <model_type type>
+void wrapIsotropicHardening(py::module& mod) {
+  constexpr UInt dim = model_type_traits<type>::dimension;
+  py::class_<IsotropicHardening<type>>(mod, "IsotropicHardening")
+      .def(py::init<Model*, Real, Real>())
+      .def_property("h", &IsotropicHardening<type>::getHardeningModulus,
+                    &IsotropicHardening<type>::setHardeningModulus)
+      .def_property("sigma_0", &IsotropicHardening<type>::getYieldStress,
+                    &IsotropicHardening<type>::setYieldStress)
+      .def("computeStress",
+           [](IsotropicHardening<type>& iso, Grid<Real, dim>& stress,
+              const Grid<Real, dim>& strain,
+              const Grid<Real, dim>& strain_increment) {
+             iso.template computeStress<false>(stress, strain, strain_increment);
+           })
+      .def("computeStressUpdate",
+           [](IsotropicHardening<type>& iso, Grid<Real, dim>& stress,
+              const Grid<Real, dim>& strain,
+              const Grid<Real, dim>& strain_increment) {
+             iso.template computeStress<true>(stress, strain, strain_increment);
+           });
+}
+
 /// Wrap temporary features for testing
 void wrapTestFeatures(py::module& mod) {
   auto test_module = mod.def_submodule("_test_features");
   test_module.doc() =
       "Module for testing new features.\n"
       "DISCLAIMER: this API is subject to frequent and unannounced changes "
       "and should **not** be relied upon!";
   wrapKOp<Kelvin, model_type::volume_2d, 2>("Kelvin_", test_module);
   wrapKOp<Kelvin, model_type::volume_2d, 3>("Kelvin_", test_module);
   wrapKOp<Kelvin, model_type::volume_2d, 4>("Kelvin_", test_module);
   wrapKOp<Mindlin, model_type::volume_2d, 3>("Mindlin_", test_module);
   wrapKOp<Mindlin, model_type::volume_2d, 4>("Mindlin_", test_module);
+
+  wrapIsotropicHardening<model_type::volume_2d>(test_module);
 }
 }
 
 __END_TAMAAS__
diff --git a/src/SConscript b/src/SConscript
index 82b287c..3256916 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,115 +1,116 @@
 import os
 Import('main_env')
 
 
 def prepend(path, list):
     return [os.path.join(path, x) for x in list]
 
 
 env = main_env
 
 # Core
 core_list = """
 fft_plan_manager.cpp
 fftransform.cpp
 fftransform_fftw.cpp
 grid.cpp
 grid_hermitian.cpp
 statistics.cpp
 surface.cpp
 tamaas.cpp
 legacy_types.cpp
 loop.cpp
 """.split()
 core_list = prepend('core', core_list)
 core_list.append('tamaas_info.cpp')
 
 # Lib roughcontact
 generator_list = """
 surface_generator.cpp
 surface_generator_filter.cpp
 surface_generator_filter_fft.cpp
 isopowerlaw.cpp
 """.split()
 generator_list = prepend('surface', generator_list)
 
 # Lib PERCOLATION
 percolation_list = """
 flood_fill.cpp
 """.split()
 percolation_list = prepend('percolation', percolation_list)
 
 # BEM PERCOLATION
 bem_list = """
 bem_kato.cpp
 bem_polonski.cpp
 bem_gigi.cpp
 bem_gigipol.cpp
 bem_penalty.cpp
 bem_uzawa.cpp
 bem_fft_base.cpp
 bem_functional.cpp
 bem_meta_functional.cpp
 elastic_energy_functional.cpp
 exponential_adhesion_functional.cpp
 squared_exponential_adhesion_functional.cpp
 maugis_adhesion_functional.cpp
 complimentary_term_functional.cpp
 bem_grid.cpp
 bem_grid_polonski.cpp
 bem_grid_kato.cpp
 bem_grid_teboulle.cpp
 bem_grid_condat.cpp
 """.split()
 bem_list = prepend('bem', bem_list)
 
 # Model
 model_list = """
 model.cpp
 model_factory.cpp
 model_type.cpp
 model_template.cpp
 be_engine.cpp
 westergaard.cpp
 elastic_functional.cpp
 meta_functional.cpp
 adhesion_functional.cpp
 elasto_plastic_model.cpp
 volume_potential.cpp
 kelvin.cpp
 mindlin.cpp
+isotropic_hardening.cpp
 """.split()
 model_list = prepend('model', model_list)
 
 # Solvers
 solvers_list = """
 contact_solver.cpp
 polonsky_keer_rey.cpp
 """.split()
 solvers_list = prepend('solvers', solvers_list)
 
 # GPU API
 gpu_list = """
 fftransform_cufft.cpp
 """.split()
 gpu_list = prepend('gpu', gpu_list)
 
 # Assembling total list
 rough_contact_list = \
   core_list + generator_list + percolation_list + model_list + solvers_list + bem_list
 
 # Adding GPU if needed
 if env['backend'] == 'cuda':
     rough_contact_list += gpu_list
 
 # Generating dependency files
 # env.AppendUnique(CXXFLAGS=['-MMD'])
 
 # for file_name in rough_contact_list:
 #     obj = file_name.replace('.cpp', '.os')
 #     dep_file_name = file_name.replace('.cpp', '.d')
 #     env.SideEffect(dep_file_name, obj)
 #     env.ParseDepends(dep_file_name)
 
 libTamaas = env.SharedLibrary('Tamaas', rough_contact_list)
 Export('libTamaas')
diff --git a/src/model/influence.hh b/src/model/influence.hh
index 7d9f217..87d5024 100644
--- a/src/model/influence.hh
+++ b/src/model/influence.hh
@@ -1,485 +1,505 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef INFLUENCE_HH
 #define INFLUENCE_HH
 /* -------------------------------------------------------------------------- */
 #include "loop.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 namespace influence {
 
 namespace {
 const Complex I(0, 1);  /// < imaginary unit
 constexpr inline Real sign(bool upper) { return (upper) ? 1. : -1.; }
 }
 
+template <UInt dim>
+struct ElasticHelper {
+  ElasticHelper(Real mu, Real nu)
+      : mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
+
+  template <typename DT, typename ST>
+  Matrix<DT, dim, dim>
+  operator()(const StaticMatrix<DT, ST, dim, dim>& gradu) const {
+    auto trace{gradu.trace()};
+    Matrix<DT, dim, dim> sigma;
+    for (UInt i = 0; i < dim; ++i)
+      for (UInt j = 0; j < dim; ++j)
+        sigma(i, j) =
+            (i == j) * lambda * trace + mu * (gradu(i, j) + gradu(j, i));
+    return sigma;
+  }
+
+  const Real mu, nu, lambda;
+};
+
 template <bool conjugate, UInt dim_q>  // dim + 1 to help dim template deduction
 inline Vector<Complex, dim_q + 1> computeD(VectorProxy<const Real, dim_q>& q) {
   constexpr Real s = -sign(conjugate);
   const auto q_norm = q.l2norm();
   return {{{s * I * q(0) / q_norm, s * I * q(1) / q_norm, 1}}};
 }
 
 template <UInt dim_q>  // dim + 1 to help dim template deduction
 inline Vector<Complex, dim_q + 1> computeD2(VectorProxy<const Real, dim_q>& q) {
   const auto q_norm = q.l2norm();
   return {{{I * q(0) / q_norm, I * q(1) / q_norm, 0}}};
 }
 
 /// Class for the Kelvin tensor
 template <UInt dim, UInt derivative_order>
 class Kelvin;
 
 /// See details in notebook and Marc Bonnet's notes
 template <>
 class Kelvin<3, 0> {
 protected:
   constexpr static UInt dim = 3;
 
 public:
   Kelvin(Real mu, Real nu) : mu(mu), b(4 * (1 - nu)) {}
 
   template <bool upper, typename ST>
   inline Vector<Complex, dim> applyU0(VectorProxy<const Real, dim - 1>& q,
                                       StaticVector<Complex, ST, dim>& f) const;
 
   template <bool upper, typename ST>
   inline Vector<Complex, dim> applyU1(VectorProxy<const Real, dim - 1>& q,
                                       StaticVector<Complex, ST, dim>& f) const;
 
 protected:
   const Real mu, b;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Base Kelvin tensor */
 /* -------------------------------------------------------------------------- */
 
 template <bool upper, typename ST>  //< this is ignored
 inline Vector<Complex, Kelvin<3, 0>::dim>
 Kelvin<3, 0>::applyU0(VectorProxy<const Real, dim - 1>& q,
                       StaticVector<Complex, ST, dim>& f) const {
   auto d2{computeD2(q)};
   auto tmp{d2.dot(f) * d2};
   tmp += b * f;
   tmp(2) -= f(2);
   tmp *= 1. / (2. * mu * q.l2norm() * b);
   return tmp;
 }
 
 template <bool upper, typename ST>
 inline Vector<Complex, Kelvin<3, 0>::dim>
 Kelvin<3, 0>::applyU1(VectorProxy<const Real, dim - 1>& q,
                       StaticVector<Complex, ST, dim>& f) const {
   auto d{computeD<upper>(q)};
   auto tmp{d.dot(f) * d};
   tmp *= sign(upper) / (2. * mu * q.l2norm() * b);
   return tmp;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Kelvin first derivative */
 /* -------------------------------------------------------------------------- */
 
 template <>
 class Kelvin<3, 1> : protected Kelvin<3, 0> {
   using parent = Kelvin<3, 0>;
 
 protected:
   static constexpr UInt dim = parent::dim;
 
 public:
   using parent::parent;
   template <bool upper, typename ST>
   inline Vector<Complex, dim>
   applyU0(VectorProxy<const Real, dim - 1>& q,
           StaticMatrix<Complex, ST, dim, dim>& f) const;
   template <bool upper, typename ST>
   inline Vector<Complex, dim>
   applyU1(VectorProxy<const Real, dim - 1>& q,
           StaticMatrix<Complex, ST, dim, dim>& f) const;
 };
 
 template <bool upper, typename ST>
 inline Vector<Complex, Kelvin<3, 1>::dim>
 Kelvin<3, 1>::applyU0(VectorProxy<const Real, dim - 1>& q,
                       StaticMatrix<Complex, ST, dim, dim>& f) const {
   auto tmp{f * computeD<upper>(q)};
   auto res{Kelvin<3, 0>::applyU0<false>(q, tmp)};
   res *= -sign(upper);
   Vector<Real, dim> e3{{{0, 0, 1}}};
   tmp.template mul<false>(f, e3);
   res += Kelvin<3, 0>::applyU1<upper>(q, tmp);
   res *= q.l2norm();
   return res;
 }
 
 template <bool upper, typename ST>
 inline Vector<Complex, Kelvin<3, 1>::dim>
 Kelvin<3, 1>::applyU1(VectorProxy<const Real, dim - 1>& q,
                       StaticMatrix<Complex, ST, dim, dim>& f) const {
   auto tmp{f * computeD<upper>(q)};
   auto res{Kelvin<3, 0>::applyU1<upper>(q, tmp)};
   res *= -sign(upper) * q.l2norm();
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Kelvin second derivative */
 /* -------------------------------------------------------------------------- */
 template <>
 class Kelvin<3, 2> : protected Kelvin<3, 1> {
   using parent = Kelvin<3, 1>;
   static constexpr UInt dim = parent::dim;
 
 public:
   using parent::parent;
   template <bool upper, typename ST>
   inline Matrix<Complex, dim, dim>
   applyU0(VectorProxy<const Real, dim - 1>& q,
           StaticMatrix<Complex, ST, dim, dim>& f) const;
   template <bool upper, typename ST>
   inline Matrix<Complex, dim, dim>
   applyU1(VectorProxy<const Real, dim - 1>& q,
           StaticMatrix<Complex, ST, dim, dim>& f) const;
 
   template <typename ST>
   inline Matrix<Complex, dim, dim>
   applyDiscontinuityTerm(VectorProxy<const Real, dim - 1>& q,
                          StaticMatrix<Complex, ST, dim, dim>& f) const;
 };
 
 template <bool upper, typename ST>
 inline Matrix<Complex, Kelvin<3, 2>::dim, Kelvin<3, 2>::dim>
 Kelvin<3, 2>::applyU0(VectorProxy<const Real, dim - 1>& q,
                       StaticMatrix<Complex, ST, dim, dim>& f) const {
   auto tmp{Kelvin<3, 1>::applyU0<upper>(q, f)};
   tmp *= -sign(upper);
   auto res{outer(computeD<upper>(q), tmp)};
   Vector<Real, dim> e3{{{0, 0, 1}}};
   res += outer(e3, Kelvin<3, 1>::applyU1<upper>(q, f));
   res *= -1. * q.l2norm();  // << -1 for grad_x
   return res;
 }
 
 template <bool upper, typename ST>
 inline Matrix<Complex, Kelvin<3, 2>::dim, Kelvin<3, 2>::dim>
 Kelvin<3, 2>::applyU1(VectorProxy<const Real, dim - 1>& q,
                       StaticMatrix<Complex, ST, dim, dim>& f) const {
   auto tmp{Kelvin<3, 1>::applyU1<upper>(q, f)};
   auto res{outer(computeD<upper>(q), tmp)};
   res *= -1. * -sign(upper) * q.l2norm();  // << -1 for grad_x
   return res;
 }
 
 template <typename ST>
 inline Matrix<Complex, Kelvin<3, 2>::dim, Kelvin<3, 2>::dim>
 Kelvin<3, 2>::applyDiscontinuityTerm(
     VectorProxy<const Real, dim - 1>& q,
     StaticMatrix<Complex, ST, dim, dim>& f) const {
   Vector<Real, dim> e3{{{0, 0, 1}}};
   const auto d2{computeD2(q)};
   auto tmp{f * e3};
   auto res{Kelvin<3, 0>::applyU0<false>(q, tmp)};
   res *= -2 * q.l2norm();
   auto scal{tmp.dot(d2)};
   auto last_component = tmp(2);
   tmp = d2;
   tmp *= scal;
   tmp(2) += last_component;
   tmp *= 1. / (this->mu * this->b);
   res += tmp;
   e3(2) *= -1;  // << -1 for grad_x
   return outer(e3, res);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Implementation */
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt order>
 struct KelvinIntegrator {
 
   static inline Real int_exp(Real a, Real z) { return std::exp(a * z) / a; }
 
   static inline Real int_exp_z(Real a, Real z) {
     const Real a_inv = 1. / a;
     return std::exp(a * z) * (z - a_inv) * a_inv;
   }
 
   static inline Real int_exp_zz(Real a, Real z) {
     const Real a_inv = 1. / a;
     return std::exp(a * z) * (z * z - 2 * a_inv * z + 2 * a_inv * a_inv) *
            a_inv;
   }
 
   template <bool upper>
   static inline Real g0(Real z) {
     return std::exp(-sign(upper) * z);
   }
 
   template <bool upper>
   static inline Real int_g0(const thrust::pair<Real, Real>& bounds,
                             Real alpha) {
     alpha *= -sign(upper);
     return int_exp(alpha, bounds.second) - int_exp(alpha, bounds.first);
   }
 
   template <bool upper>
   static inline Real int_g0_z(const thrust::pair<Real, Real>& bounds,
                               Real alpha) {
     alpha *= -sign(upper);
     return int_exp_z(alpha, bounds.second) - int_exp_z(alpha, bounds.first);
   }
 
   template <bool upper>
   static inline Real g1(Real z) {
     return g0<upper>(z) * z;
   }
 
   template <bool upper>
   static inline Real int_g1(const thrust::pair<Real, Real>& bounds,
                             Real alpha) {
     return alpha * int_g0_z<upper>(bounds, alpha);
   }
 
   template <bool upper>
   static inline Real int_g1_z(const thrust::pair<Real, Real>& bounds,
                               Real alpha) {
     return alpha * (int_exp_zz(-alpha * sign(upper), bounds.second) -
                     int_exp_zz(-alpha * sign(upper), bounds.first));
   }
 
   template <Int yj_xi, UInt dim_q, typename OutType, typename InType,
             typename KelvinType>
   static inline void integrate(OutType&& out, InType&& in, KelvinType&& kelvin,
                                VectorProxy<const Real, dim_q>& q_, Real dij,
                                Real dl);
 };
 
 // Integrate against a constant shape function
 template <>
 template <Int yj_xi, UInt dim_q, typename OutType, typename InType,
           typename KelvinType>
 inline void KelvinIntegrator<0>::integrate(OutType&& out, InType&& in,
                                            KelvinType&& kelvin,
                                            VectorProxy<const Real, dim_q>& q_,
                                            Real dij, Real dl) {
   const auto q = q_.l2norm();
   const auto alpha = q * dl, beta = q * dij;
 
   // Integrate on [-1, 0]
   const auto bounds_low = thrust::make_pair(-0.5, 0);
   constexpr bool domain_low = yj_xi > 0;
 
   auto g0_mul{kelvin.template applyU0<domain_low>(q_, in)};
   auto g1_mul{kelvin.template applyU1<domain_low>(q_, in)};
 
   decltype(g0_mul) res;
 
   res = g0_mul * (g0<domain_low>(beta) * int_g0<domain_low>(bounds_low, alpha));
   res +=
       g1_mul * (g1<domain_low>(beta) * int_g0<domain_low>(bounds_low, alpha) +
                 g0<domain_low>(beta) * int_g1<domain_low>(bounds_low, alpha));
 
   // Integrate on [0, 1]
   const auto bounds_up = thrust::make_pair(0, 0.5);
   constexpr bool domain_up = yj_xi >= 0;
 
   // Avoid computing
   if (yj_xi == 0) {
     g0_mul = kelvin.template applyU0<true>(q_, in);
     g1_mul = kelvin.template applyU1<true>(q_, in);
   }
 
   res += g0_mul * (g0<domain_up>(beta) * int_g0<domain_up>(bounds_up, alpha));
   res += g1_mul * (g1<domain_up>(beta) * int_g0<domain_up>(bounds_up, alpha) +
                    g0<domain_up>(beta) * int_g1<domain_up>(bounds_up, alpha));
 
   res *= dl;
   out += res;
 }
 
 // Integrate against a linear shape function
 template <>
 template <Int yj_xi, UInt dim_q, typename OutType, typename InType,
           typename KelvinType>
 inline void KelvinIntegrator<1>::integrate(OutType&& out, InType&& in,
                                            KelvinType&& kelvin,
                                            VectorProxy<const Real, dim_q>& q_,
                                            Real dij, Real dl) {
   const auto q = q_.l2norm();
   const auto alpha = q * dl, beta = q * dij;
 
   // Integrate on [-1, 0]
   const auto bounds_low = thrust::make_pair(-1, 0);
   constexpr bool domain_low = yj_xi > 0;
 
   auto g0_mul{kelvin.template applyU0<domain_low>(q_, in)};
   auto g1_mul{kelvin.template applyU1<domain_low>(q_, in)};
 
   decltype(g0_mul) res;
 
   res = g0_mul *
         (g0<domain_low>(beta) * (int_g0<domain_low>(bounds_low, alpha) +
                                  int_g0_z<domain_low>(bounds_low, alpha)));
   res += g1_mul *
          (g1<domain_low>(beta) * (int_g0<domain_low>(bounds_low, alpha) +
                                   int_g0_z<domain_low>(bounds_low, alpha)) +
           g0<domain_low>(beta) * (int_g1<domain_low>(bounds_low, alpha) +
                                   int_g1_z<domain_low>(bounds_low, alpha)));
 
   // Integrate on [0, 1]
   const auto bounds_up = thrust::make_pair(0, 1);
   constexpr bool domain_up = yj_xi >= 0;
 
   // Avoid computing
   if (yj_xi == 0) {
     g0_mul = kelvin.template applyU0<true>(q_, in);
     g1_mul = kelvin.template applyU1<true>(q_, in);
   }
 
   res +=
       g0_mul * (g0<domain_up>(beta) * (int_g0<domain_up>(bounds_up, alpha) -
                                        int_g0_z<domain_up>(bounds_up, alpha)));
   res +=
       g1_mul * (g1<domain_up>(beta) * (int_g0<domain_up>(bounds_up, alpha) -
                                        int_g0_z<domain_up>(bounds_up, alpha)) +
                 g0<domain_up>(beta) * (int_g1<domain_up>(bounds_up, alpha) -
                                        int_g1_z<domain_up>(bounds_up, alpha)));
 
   res *= dl;
   out += res;
 }
 
 /* -------------------------------------------------------------------------- */
 /// Class for the Boussinesq tensor
 template <UInt dim, UInt derivative_order>
 class Boussinesq;
 
 template <>
 class Boussinesq<3, 0> {
 protected:
   static constexpr UInt dim = 3;
 
 public:
   /// Constructor
   Boussinesq(Real mu, Real nu)
       : mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
 
   template <typename ST>
   inline Vector<Complex, dim>
   applyU0(StaticVector<Complex, ST, dim>& t,
           VectorProxy<const Real, dim - 1>& q) const;
   template <typename ST>
   inline Vector<Complex, dim>
   applyU1(StaticVector<Complex, ST, dim>& t,
           VectorProxy<const Real, dim - 1>& q) const;
 
 protected:
   const Real mu, nu, lambda;
 };
 
 template <typename ST>
 inline Vector<Complex, Boussinesq<3, 0>::dim>
 Boussinesq<3, 0>::applyU0(StaticVector<Complex, ST, dim>& t,
                           VectorProxy<const Real, dim - 1>& q) const {
   // anticipating K(-q) convolution (no change for d2 cause d2xd2 not sensitive)
   auto d{computeD<true>(q)}, dbar{computeD<false>(q)}, d2{computeD2(q)};
   auto res{2. * t};
   res += ((1. - 2. * nu) * dbar.dot(t)) * d;
   res += d2.dot(t) * d2;
   res(2) -= t(2);
   res *= (1. / (2. * mu * q.l2norm()));
   return res;
 }
 
 template <typename ST>
 inline Vector<Complex, Boussinesq<3, 0>::dim>
 Boussinesq<3, 0>::applyU1(StaticVector<Complex, ST, dim>& t,
                           VectorProxy<const Real, dim - 1>& q) const {
   // anticipating K(-q) convolution (no change for d2 cause d2xd2 not sensitive)
   auto dbar{computeD<false>(q)};
   auto res{(dbar.dot(t) / (2. * mu * q.l2norm())) * dbar};
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 class Boussinesq<3, 1> : protected Boussinesq<3, 0> {
 protected:
   using parent = Boussinesq<3, 0>;
   static const UInt dim = parent::dim;
 
 public:
   using parent::parent;
   template <typename ST>
   inline Matrix<Complex, dim, dim>
   applyU0(StaticVector<Complex, ST, dim>& t,
           VectorProxy<const Real, dim - 1>& q) const;
 
   template <typename ST>
   inline Matrix<Complex, dim, dim>
   applyU1(StaticVector<Complex, ST, dim>& t,
           VectorProxy<const Real, dim - 1>& q) const;
 };
 
 template <typename ST>
 inline Matrix<Complex, Boussinesq<3, 1>::dim, Boussinesq<3, 1>::dim>
 Boussinesq<3, 1>::applyU0(StaticVector<Complex, ST, dim>& t,
                           VectorProxy<const Real, dim - 1>& q) const {
   // taking into account K(-q) convolution
   auto dbar{computeD<false>(q)};
   Vector<Real, dim> e3{{{0, 0, 1}}};
   auto res{outer(dbar, Boussinesq<3, 0>::applyU0(t, q))};
   res *= -1.;
   res += outer(e3, Boussinesq<3, 0>::applyU1(t, q));
   res *= q.l2norm();
   return res;
 }
 
 template <typename ST>
 inline Matrix<Complex, Boussinesq<3, 1>::dim, Boussinesq<3, 1>::dim>
 Boussinesq<3, 1>::applyU1(StaticVector<Complex, ST, dim>& t,
                           VectorProxy<const Real, dim - 1>& q) const {
   // taking into account K(-q) convolution
   auto dbar{computeD<false>(q)};
   auto res{outer(dbar, Boussinesq<3, 0>::applyU1(t, q))};
   res *= -q.l2norm();
   return res;
 }
 }  // namespace influence
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 
 #endif  // INFLUENCE_HH
diff --git a/src/model/isotropic_hardening.cpp b/src/model/isotropic_hardening.cpp
new file mode 100644
index 0000000..9ddd24b
--- /dev/null
+++ b/src/model/isotropic_hardening.cpp
@@ -0,0 +1,48 @@
+/**
+ * @file
+ *
+ * @author Lucas Frérot <lucas.frerot@epfl.ch>
+ *
+ * @section LICENSE
+ *
+ * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
+ * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
+ * Solides)
+ *
+ * Tamaas is free  software: you can redistribute it and/or  modify it under the
+ * terms  of the  GNU Lesser  General Public  License as  published by  the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Tamaas 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  Lesser General  Public License  for more
+ * details.
+ *
+ * You should  have received  a copy  of the GNU  Lesser General  Public License
+ * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#include "isotropic_hardening.hh"
+#include "influence.hh"
+/* -------------------------------------------------------------------------- */
+__BEGIN_TAMAAS__
+/* -------------------------------------------------------------------------- */
+
+template <model_type type>
+IsotropicHardening<type>::IsotropicHardening(Model* model, Real sigma_0, Real h)
+    : model(model), sigma_0(sigma_0), h(h) {
+  plastic_strain.setNbComponents(trait::components * trait::components);
+  cumulated_plastic_strain.setNbComponents(1);
+
+  plastic_strain.resize(model->getDiscretization());
+  cumulated_plastic_strain.resize(model->getDiscretization());
+}
+
+/* -------------------------------------------------------------------------- */
+/* Template instanciation */
+/* -------------------------------------------------------------------------- */
+template class IsotropicHardening<model_type::volume_2d>;
+
+__END_TAMAAS__
diff --git a/src/model/isotropic_hardening.hh b/src/model/isotropic_hardening.hh
new file mode 100644
index 0000000..66fae9b
--- /dev/null
+++ b/src/model/isotropic_hardening.hh
@@ -0,0 +1,113 @@
+/**
+ * @file
+ *
+ * @author Lucas Frérot <lucas.frerot@epfl.ch>
+ *
+ * @section LICENSE
+ *
+ * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
+ * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
+ * Solides)
+ *
+ * Tamaas is free  software: you can redistribute it and/or  modify it under the
+ * terms  of the  GNU Lesser  General Public  License as  published by  the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Tamaas 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  Lesser General  Public License  for more
+ * details.
+ *
+ * You should  have received  a copy  of the GNU  Lesser General  Public License
+ * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#ifndef ISOTROPIC_HARDENING_HH
+#define ISOTROPIC_HARDENING_HH
+/* -------------------------------------------------------------------------- */
+#include "grid.hh"
+#include "influence.hh"
+#include "model.hh"
+#include "model_type.hh"
+/* -------------------------------------------------------------------------- */
+__BEGIN_TAMAAS__
+/* -------------------------------------------------------------------------- */
+
+template <model_type type>
+class IsotropicHardening {
+  using trait = model_type_traits<type>;
+  static constexpr UInt dim = trait::dimension;
+
+public:
+  /// Constructor
+  IsotropicHardening(Model* model, Real sigma_0, Real h);
+
+  /// Compute stresses with radial return algorithm
+  template <bool update>
+  void computeStress(Grid<Real, dim>& stress, const Grid<Real, dim>& strain,
+                     const Grid<Real, dim>& strain_increment);
+
+  Real hardening(Real p) { return sigma_0 + h * p; }
+
+  TAMAAS_ACCESSOR(h, Real, HardeningModulus);
+  TAMAAS_ACCESSOR(sigma_0, Real, YieldStress);
+
+protected:
+  Model* model;
+  Real sigma_0;  /// < initial yield stress
+  Real h;        /// < hardening modulus
+  Grid<Real, dim> plastic_strain, cumulated_plastic_strain;
+};
+
+/* -------------------------------------------------------------------------- */
+template <>
+template <bool update>
+void IsotropicHardening<model_type::volume_2d>::computeStress(
+    Grid<Real, dim>& stress, const Grid<Real, dim>& strain,
+    const Grid<Real, dim>& strain_increment) {
+  influence::ElasticHelper<dim> elasticity(this->model->getShearModulus(),
+                                           this->model->getPoissonRatio());
+  Loop::stridedLoop(
+      [&elasticity,
+       this](MatrixProxy<Real, dim, dim>&& sigma,          // < projected stress
+             MatrixProxy<const Real, dim, dim>&& epsilon,  // < total strain
+             MatrixProxy<const Real, dim, dim>&&
+                 delta_epsilon,                 // < strain increment
+             MatrixProxy<Real, dim, dim>&& ep,  // < plastic strain
+             Real& p /* < cumulated plastic strain */) {
+        auto sigma_tr{elasticity(epsilon + delta_epsilon)};
+        sigma = sigma_tr;
+
+        Matrix<Real, dim, dim> dev;
+        const Real hydrostatic_stress{sigma_tr.trace() / 3.};
+
+        for (UInt i = 0; i < dim; ++i)
+          for (UInt j = 0; j < dim; ++j)
+            dev(i, j) = sigma_tr(i, j) - (i == j) * hydrostatic_stress;
+
+        Real von_mises{std::sqrt(1.5) * dev.l2norm()};
+
+        if (von_mises - hardening(p) > 0) {
+          // radial return
+          Real dp = (von_mises - hardening(p)) / (3 * elasticity.mu + h);
+
+          dev *= 3 * dp / (2 * von_mises);  // saving memory (dev is delta ep)
+
+          if (update) {  // constexpr if when C++17
+            p += dp;
+            ep += dev;
+          }
+
+          dev *= 2 * elasticity.mu;
+          sigma -= dev;
+        }
+      },
+      stress, strain, strain_increment, this->plastic_strain,
+      this->cumulated_plastic_strain);
+}
+
+__END_TAMAAS__
+/* -------------------------------------------------------------------------- */
+#endif  // ISOTROPIC_HARDENING_HH
diff --git a/src/model/mindlin.cpp b/src/model/mindlin.cpp
index 7d0bbf5..f159fd2 100644
--- a/src/model/mindlin.cpp
+++ b/src/model/mindlin.cpp
@@ -1,186 +1,161 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "mindlin.hh"
 #include "elasto_plastic_model.hh"
 #include "influence.hh"
 #include "kelvin_helper.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
-namespace influence {
-template <UInt dim>
-class ElasticHelper {
-public:
-  ElasticHelper(Real mu, Real nu)
-      : mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
-
-  template <typename ST>
-  Matrix<Complex, dim, dim>
-  operator()(StaticMatrix<Complex, ST, dim, dim>& gradu) const {
-    auto trace{gradu.trace()};
-    Matrix<Complex, dim, dim> sigma;
-    for (UInt i = 0; i < dim; ++i)
-      for (UInt j = 0; j < dim; ++j)
-        sigma(i, j) =
-            (i == j) * lambda * trace + mu * (gradu(i, j) + gradu(j, i));
-    return sigma;
-  }
-
-protected:
-  const Real mu, nu, lambda;
-};
-
-}  // namespace influence
-
 namespace detail {
 template <model_type type, typename T>
 class MindlinBoussinesqHelper {
   using boussinesq_t = T;
   using trait = model_type_traits<type>;
   static constexpr UInt dim = trait::dimension;
   static constexpr UInt bdim = trait::boundary_dimension;
 
 public:
   MindlinBoussinesqHelper(const T& boussinesq,
                           const influence::ElasticHelper<dim>& el, Real y_3,
                           Real cutoff)
       : boussinesq(boussinesq), el(el), y_3(y_3), cutoff(cutoff) {}
 
   inline void applyIntegral(GridHermitian<Real, bdim>& out,
                             GridHermitian<Real, bdim>& source,
                             const Grid<Real, bdim>& wavevectors) const {
     Loop::stridedLoop(
         [&](typename KelvinTrait<T>::out_t&& out_local,
             MatrixProxy<Complex, dim, dim>&& source_local,
             VectorProxy<const Real, bdim>&& q) {
           // Cutoff
           if (-q.l2norm() * std::abs(y_3) < std::log(cutoff))
             return;
 
 	  Vector<Real, dim> normal{{{0, 0, 1}}};
 	  auto traction{el(source_local) * normal};
 	  out_local += boussinesq.applyU0(traction, q) *
 	    influence::KelvinIntegrator<0>::g0<true>(q.l2norm() * y_3);
 	  out_local += (boussinesq.applyU1(traction, q) *
 		  influence::KelvinIntegrator<0>::g1<true>(q.l2norm() * y_3));
         },
         out, source, wavevectors);
   }
 
 protected:
   const T& boussinesq;
   const influence::ElasticHelper<dim>& el;
   const Real y_3, cutoff;
 };
 }  // namespace detail
 
 template <model_type type, UInt order>
 void Mindlin<type, order>::apply(GridBase<Real>& source,
                                  GridBase<Real>& out) const {
   Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus();
   constexpr UInt derivative = order - 2;
   influence::Kelvin<trait::dimension, derivative> kelvin(mu, nu);
   influence::Kelvin<trait::dimension, 2> kelvin_strain(mu, nu);
   influence::Boussinesq<trait::dimension, derivative - 1> boussinesq(mu, nu);
   influence::ElasticHelper<trait::dimension> elasticity(mu, nu);
 
   auto apply = [&](UInt i, decltype(this->source_buffers)& source_buffers,
                    decltype(this->disp_buffer)& out_buffer) {
     constexpr UInt dim = trait::dimension;
     const Real L = this->model->getSystemSize().front();
     const UInt N = this->model->getDiscretization().front();
     const Real dl = L / (N - 1);
 
     out_buffer = 0;
 
     /// Compute surface strains only once (dirty)
     if (i == 0) {
       surface_strains = 0;
       for (UInt j : Loop::range(N)) {
         const Real dij = j * dl;
         auto& source = source_buffers[j];
         detail::KelvinHelper<type, influence::Kelvin<dim, 2>> helper(
             kelvin_strain, dij, dl, 1e-2);
 
         if (j > i) {
           helper.template applyIntegral<1>(surface_strains, source,
                                            this->wavevectors);
         } else if (j == i) {
           helper.template applyIntegral<0>(surface_strains, source,
                                            this->wavevectors);
           helper.applyFreeTerm(surface_strains, source, this->wavevectors);
         } else {
           helper.template applyIntegral<-1>(surface_strains, source,
                                             this->wavevectors);
         }
       }
     }
 
     for (UInt j : Loop::range(N)) {
       const Real dij = j * dl - i * dl;  // don't factorize!
       auto& source = source_buffers[j];
       detail::KelvinHelper<type, influence::Kelvin<dim, derivative>> helper(
           kelvin, dij, dl, 1e-2);
 
       if (j > i) {
         helper.template applyIntegral<1>(out_buffer, source, this->wavevectors);
       } else if (j == i) {
         helper.template applyIntegral<0>(out_buffer, source, this->wavevectors);
         helper.applyFreeTerm(out_buffer, source, this->wavevectors);
       } else {
         helper.template applyIntegral<-1>(out_buffer, source,
                                           this->wavevectors);
       }
     }
 
     // Correcting for the tractions on the surface
     Real xi = i * dl;
     detail::MindlinBoussinesqHelper<type,
                                     influence::Boussinesq<dim, derivative - 1>>
         helper(boussinesq, elasticity, xi, 1e-2);
     helper.applyIntegral(out_buffer, surface_strains, this->wavevectors);
     
 
     // Setting fundamental frequency to zero
     typename detail::KelvinTrait<
         influence::Kelvin<trait::dimension, derivative>>::out_t
         out_fundamental(out_buffer(0));
     out_fundamental = 0;
   };
 
   this->fourierApply(apply, source, out);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation                                                     */
 /* -------------------------------------------------------------------------- */
 template class Mindlin<model_type::volume_2d, 3>;
 template class Mindlin<model_type::volume_2d, 4>;
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__