diff --git a/python/tamaas/nonlinear_solvers/__init__.py b/python/tamaas/nonlinear_solvers/__init__.py
index b96b06c..a23627c 100644
--- a/python/tamaas/nonlinear_solvers/__init__.py
+++ b/python/tamaas/nonlinear_solvers/__init__.py
@@ -1,221 +1,221 @@
 # -*- mode:python; coding: utf-8 -*-
 # @file
 # @section LICENSE
 #
 # Copyright (©) 2016-2020 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/>.
 
 """
 Pulling solvers to nonlinear_solvers module
 """
 
 from .. import EPSolver
 
 import numpy as np
 from scipy.sparse.linalg import LinearOperator, lgmres
 from scipy.linalg import norm
 from scipy.optimize import newton_krylov, root
 
 try:
     from .akantu_solver import *  # noqa
 except ImportError:
     pass
 
 __all__ = ['NLNoConvergence',
            'DFSANESolver',
            'NewtonKrylovSolver',
            'ToleranceManager']
 
 
 class NLNoConvergence(Exception):
     """Convergence not reached exception"""
     pass
 
 
 class NewtonRaphsonSolver(EPSolver):
     """
     Newton-Raphson based on CTO
     """
 
     def __init__(self, residual, model):
         super(NewtonRaphsonSolver, self).__init__(residual)
         self.model = model
 
         self._residual = self.getResidual()
         self._x = self.getStrainIncrement()
         self.last_solution = np.zeros_like(self._x)
 
     def applyTangent(self, _input, strain):
         _input = np.reshape(_input, strain.shape)
         out = np.zeros_like(_input)
         self._residual.applyTangent(out, _input, strain)
         return out
 
     def solve(self):
         self._x[...] = 0
 
         # For initial guess, compute the strain due to boundary tractions
         self._residual.computeResidual(self._x)
         self._x[...] = self._residual.getVector()
 
         res_vec = np.ravel(self._residual.getVector())
         initial_norm = norm(res_vec)
         self._residual.computeResidual(self._x)
 
         n_max = 100
         i = 0
 
         while norm(res_vec) / initial_norm > 1e-10 and i < n_max:
             # Making linear tangent
             tangent = LinearOperator(
                 (self._x.size, self._x.size),
                 matvec=lambda x: self.applyTangent(x, self._x))
             res_vec *= -1
             correction, ok = lgmres(tangent, res_vec, tol=1e-12, maxiter=100)
             if ok != 0:
                 raise Exception("LGMRES not converged")
             self._x += np.reshape(correction, self._x.shape)
             self._residual.computeResidual(self._x)
             i += 1
 
         print("Solved Newton-Raphson in {} iterations".format(i))
 
         # Computing the strain correction
         self.last_solution *= -1
         self.last_solution += self._x
 
         # Computing displacements
         self._residual.computeResidualDisplacement(self.last_solution)
         self.last_solution[...] = self._x[...]
 
 
 class ScipySolver(EPSolver):
     """
     Base class for solvers wrapping SciPy routines
     """
 
     def __init__(self, residual, model, callback=None):
         super(ScipySolver, self).__init__(residual)
         self.model = model
         self.callback = callback
         self._x = self.getStrainIncrement()
         self._residual = self.getResidual()
         self.options = {'ftol': 0, 'fatol': 1e-9}
 
     @property
     def tolerance(self):
         return self.options['fatol']
 
     @tolerance.setter
-    def set_tolerance(self, v):
+    def tolerance(self, v):
         self.options['fatol'] = v
 
     def solve(self):
         """
         Solve the nonlinear plasticity equation using the scipy_solve routine
         """
         # For initial guess, compute the strain due to boundary tractions
         # self._residual.computeResidual(self._x)
         # self._x[...] = self._residual.getVector()
 
         # Scipy root callback
         def compute_residual(x):
             self._residual.computeResidual(x)
             return self._residual.getVector().copy()
 
         # Solve
         self._x[...] = self.scipy_solve(compute_residual)
 
         # Computing displacements
         self._residual.computeResidualDisplacement(self._x)
 
     def reset(self):
         self._x[...] = 0
 
 
 class NewtonKrylovSolver(ScipySolver):
     """
     Solve using a finite-difference Newton-Krylov method
     """
     def __init__(self, residual, model, callback=None):
         ScipySolver.__init__(self, residual, model, callback)
 
     def scipy_solve(self, compute_residual):
         try:
             return newton_krylov(compute_residual, self._x,
                                  f_tol=self.options['fatol'],
                                  verbose=True, callback=self.callback)
         except Exception as e:
             raise NLNoConvergence(e.what())
 
 
 class DFSANESolver(ScipySolver):
     """
     Solve using a spectral residual jacobianless method
     """
     def __init__(self, residual, model, callback=None):
         ScipySolver.__init__(self, residual, model, callback)
 
     def scipy_solve(self, compute_residual):
         solution = root(compute_residual,
                         self._x,
                         method='df-sane',
                         options=self.options,
                         callback=self.callback)
         print("DF-SANE: {} ({} iterations, {})".format(
             solution.message,
             solution.nit,
             self.options))
 
         if not solution.success:
             raise NLNoConvergence("DF-SANE did not converge")
         return solution.x.copy()
 
 
 class ToleranceManager:
     """
     Decorator to manage tolerance of non-linear solver
     """
     def __init__(self, start, end, rate, key='fatol'):
         self.start = start / rate  # just anticipating first multiplication
         self.end = end
         self.rate = rate
         self.key = key
 
     def __call__(self, cls):
         orig_init = cls.__init__
         orig_solve = cls.scipy_solve
         orig_updateState = cls.updateState
 
         def __init__(obj, *args, **kwargs):
             orig_init(obj, *args, **kwargs)
             obj.options[self.key] = self.start
 
         def scipy_solve(obj, *args, **kwargs):
             ftol = obj.options[self.key]
             ftol *= self.rate
 
             obj.options[self.key] = max(ftol, self.end)
             return orig_solve(obj, *args, **kwargs)
 
         def updateState(obj, *args, **kwargs):
             obj.options[self.key] = self.start
             return orig_updateState(obj, *args, **kwargs)
 
         cls.__init__ = __init__
         cls.scipy_solve = scipy_solve
         cls.updateState = updateState
         return cls
diff --git a/tests/conftest.py b/tests/conftest.py
index d098de0..c83ff62 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,234 +1,238 @@
 # -*- coding: utf-8 -*-
 # @file
 # @section LICENSE
 #
 # Copyright (©) 2016-2020 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/>.
 
 from __future__ import division, print_function
 
 import numpy as np
 import tamaas as tm
 import pytest
 
 from numpy.linalg import norm
 
 
 def profile(func, size, mode, amplitude):
     x = np.linspace(0, 1, size[0], endpoint=False, dtype=tm.dtype)
     y = np.linspace(0, 1, size[1], endpoint=False, dtype=tm.dtype)
     x, y = np.meshgrid(x, y, indexing='ij')
     surface = amplitude * func(2*np.pi*x*mode[0]) * func(2*np.pi*y*mode[1])
     return surface.copy()
 
 
 @pytest.fixture(scope="session")
 def tamaas_fixture():
     tm.initialize()
     yield None
     tm.finalize()
 
 
 class HertzFixture:
     def __init__(self, n, load):
         self.domain_size = 1
         self.n = n
         self.load = load
         self.curvature = 0.1
         self.radius = 1. / self.curvature
         self.e_star = 1.
         self.a = (3 * load / (4 * self.curvature * self.e_star))**(1. / 3.)
 
         self.x = np.linspace(-self.domain_size / 2.,
                              self.domain_size / 2.,
                              self.n, dtype=tm.dtype)
         self.y = self.x.copy()
         self.x, self.y = np.meshgrid(self.x, self.y)
 
         self._computeSurface()
         self._computePressure()
         self._computeDisplacement()
 
     def _computeDisplacement(self):
         r = np.sqrt(self.x**2 + self.y**2)
         self.displacement = np.zeros_like(r)
         contact = r < self.a
         self.displacement[contact] = self.surface[contact]
         self.displacement[~contact] = \
             (self.surface[~contact]
              + self.a / (np.pi * self.radius) * np.sqrt(r[~contact]**2
                                                         - self.a**2)
              + (r[~contact]**2 - 2 * self.a**2) / (np.pi * self.radius)
              * np.arccos(self.a / r[~contact]))
 
     def _computePressure(self):
         r = np.sqrt(self.x**2 + self.y**2)
         self.pressure = np.zeros_like(r)
         contact = np.where(r < self.a)
         self.pressure[contact] = \
             2 * self.e_star / (np.pi * self.radius) \
             * np.sqrt(self.a**2 - r[contact]**2)
 
     def _computeSurface(self):
         self.surface = -1. / (2 * self.radius) * (self.x**2 + self.y**2)
 
 
 @pytest.fixture(scope="module")
 def hertz(tamaas_fixture):
     return HertzFixture(1024, 0.00001)
 
 
 @pytest.fixture(scope="module")
 def hertz_coarse(tamaas_fixture):
     return HertzFixture(512, 0.0001)
 
 @pytest.fixture(scope="module",
                 params=[tm.PolonskyKeerRey.pressure])
 def pkr(hertz, request):
     model = tm.ModelFactory.createModel(tm.model_type.basic_2d,
                                         [hertz.domain_size, hertz.domain_size],
                                         [hertz.n, hertz.n])
     model.setElasticity(hertz.e_star, 0)
     solver = tm.PolonskyKeerRey(model, hertz.surface, 1e-12,
                                 request.param,
                                 request.param)
     solver.solve(hertz.load)
 
     return model, hertz
 
 
 class WestergaardFixture:
     def __init__(self, n, load):
         self.domain_size = 1.
         self.lamda = 1.
         self.delta = 0.1
         self.e_star = 1.
 
         self.n = n
 
         self.p_star = np.pi * self.e_star * self.delta / self.lamda
 
         self.load = load * self.p_star
 
         self.a = self.lamda / np.pi \
             * np.arcsin(np.sqrt(self.load / self.p_star))
 
         self.x = np.linspace(-self.domain_size / 2.,
                              self.domain_size / 2.,
                              self.n, endpoint=False, dtype=tm.dtype)
 
         self._computeSurface()
         self._computePressure()
         self._computeDisplacement()
 
     def _computeSurface(self):
         self.surface = self.delta * np.cos(2 * np.pi * self.x / self.lamda)
 
     def _computePressure(self):
         self.pressure = np.zeros_like(self.surface)
         contact = np.where(np.abs(self.x) < self.a)
         self.pressure[contact] = 2 * self.load \
             * (np.cos(np.pi * self.x[contact] / self.lamda)
                / np.sin(np.pi * self.a / self.lamda)**2) \
             * np.sqrt(np.sin(np.pi * self.a / self.lamda)**2 -
                       np.sin(np.pi * self.x[contact] / self.lamda)**2)
 
     def _computeDisplacement(self):
         psi = np.pi * np.abs(self.x) / self.lamda
         psi_a = np.pi * self.a / self.lamda
 
         with np.errstate(invalid='ignore'):  # get some warnings out of the way
             self.displacement = (np.cos(2*psi) + 2 * np.sin(psi)
                                  * np.sqrt(np.sin(psi)**2 - np.sin(psi_a)**2)
                                  - 2 * np.sin(psi_a)**2 *
                                  np.log((np.sin(psi) +
                                          np.sqrt(np.sin(psi)**2
                                                  - np.sin(psi_a)**2))
                                         / np.sin(psi_a)))
 
         contact = np.where(np.abs(self.x) < self.a)
         self.displacement[contact] = np.cos(2*psi[contact])
         self.displacement *= self.load * self.lamda / (np.pi * self.e_star
                                                        * np.sin(psi_a)**2)
 
 
 @pytest.fixture(scope="module")
 def westergaard(tamaas_fixture):
     return WestergaardFixture(19683, 0.1)
 
 
 class PatchWestergaard:
     def __init__(self, model_type, domain, modes, size):
         self.E = 3.
         self.nu = 0.
         self.e_star = self.E / (1 - self.nu**2)
         self.n = size
         self.domain = domain
         self.model = tm.ModelFactory.createModel(model_type, domain,
                                                  size)
         self.model.setElasticity(self.E, self.nu)
 
         self.pressure = profile(np.cos, size, modes, 1)
         self.solution = profile(np.cos, size, modes,
                                 1 / (np.pi * self.e_star * norm(modes)))
 
 
 @pytest.fixture(scope="module", params=[tm.model_type.basic_2d])
 def patch_westergaard(tamaas_fixture, request):
     return PatchWestergaard(request.param, [1., 1.], [3, 1], [6, 6])
 
 
 class UniformPlasticity:
     def __init__(self, model_type, domain, sizes):
         self.n = sizes
         self.domain = domain
         self.model = tm.ModelFactory.createModel(model_type, domain,
                                                  sizes)
         self.E_h = 0.1
         self.sigma_y = 0.01
         self.residual = tm.ModelFactory.createResidual(self.model,
                                                        sigma_y=self.sigma_y,
                                                        hardening=self.E_h)
         self.model.E = 1.
         self.model.nu = 0.
 
     def solution(self, p):
         E, nu = self.model.E, self.model.nu
         E_h, sigma_y = self.E_h, self.sigma_y
         mu = E / (2 * (1 + nu))
 
         strain = -1 / (mu + E_h) * (p * (3 * mu + E_h) / (2 * mu)
                                     - np.sign(p) * sigma_y)
 
         dep = (2 * mu * np.abs(strain) - sigma_y) / (3 * mu + E_h)
         plastic_strain = np.sign(strain) / 2 * dep * np.array([
             -1, -1, 2, 0, 0, 0,
         ], dtype=tm.dtype)
 
         stress = 2 * mu * (np.array([
             0, 0, strain, 0, 0, 0
         ], dtype=tm.dtype) - plastic_strain)
 
         return {
             "stress": stress,
             "plastic_strain": plastic_strain,
             "cumulated_plastic_strain": dep
+        }, {
+            "stress": mu,
+            "plastic_strain": 1,
+            "cumulated_plastic_strain": 1
         }
 
 
 @pytest.fixture(scope="module", params=[tm.model_type.volume_2d])
 def patch_isotropic_plasticity(tamaas_fixture, request):
     return UniformPlasticity(request.param, [1., 1., 1.], [4, 4, 4])
diff --git a/tests/test_patch_plasticity.py b/tests/test_patch_plasticity.py
index 424b951..fe0e408 100644
--- a/tests/test_patch_plasticity.py
+++ b/tests/test_patch_plasticity.py
@@ -1,42 +1,43 @@
 # -*- coding: utf-8 -*-
 # @file
 # @section LICENSE
 #
 # Copyright (©) 2016-20 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/>.
 
 from numpy.linalg import norm
 from tamaas.nonlinear_solvers import DFSANESolver
 
 
 def test_patch_plasticity(patch_isotropic_plasticity):
     model = patch_isotropic_plasticity.model
     residual = patch_isotropic_plasticity.residual
 
     applied_pressure = 0.1
 
     solver = DFSANESolver(residual, model)
+    solver.tolerance = 1e-15
     pressure = model['traction'][..., 2]
     pressure[:] = applied_pressure
 
     solver.solve()
     solver.updateState()
 
-    solution = patch_isotropic_plasticity.solution(applied_pressure)
+    solution, normal = patch_isotropic_plasticity.solution(applied_pressure)
 
     for key in solution:
-        error = norm(model[key] - solution[key]) / applied_pressure
-        assert error < 3e-15
+        error = norm(model[key] - solution[key]) / normal[key]
+        assert error < 1e-15