diff --git a/python/SConscript b/python/SConscript
index e3bd799..0ee3a02 100644
--- a/python/SConscript
+++ b/python/SConscript
@@ -1,44 +1,50 @@
 from __future__ import print_function
 
 from os.path import abspath
 from SCons.Script import Import, Depends, Split, Copy
 
 Import('main_env')
 
 # Pybind11 wrapper
 env_pybind = main_env.Clone(SHLIBPREFIX='')
 env_pybind.Tool(pybind11)
 
 pybind_sources = Split("""
 tamaas_module.cpp
 wrap/core.cpp
 wrap/percolation.cpp
 wrap/surface.cpp
 wrap/model.cpp
 wrap/solvers.cpp
 wrap/test_features.cpp
 """)
 
 if main_env['CXX'] != 'icpc':
     pybind_sources += ["wrap/bem.cpp"]
     env_pybind.AppendUnique(CPPDEFINES="LEGACY_BEM")
 
 # Building the pybind library
 tamaas_wrap = env_pybind.SharedLibrary(
     target='tamaas/_tamaas',
     source=pybind_sources,
     LIBS=['Tamaas'],
     RPATH=[abspath('../src')]
 )
 
 # For some reason link happens too early
 Import('libTamaas')
 Depends(tamaas_wrap, libTamaas)
 
 # Copying the __init__.py file with extra python classes
 copy_env = env_pybind.Clone(
     PRINT_CMD_LINE_FUNC=main_env['gen_print']("Copying", "red", main_env))
-copy_env.Command("tamaas/__init__.py", "#/python/tamaas.py",
-                 Copy("$TARGET", "$SOURCE"))
-copy_env.Command("setup.py", "#/python/setup.py",
-                 Copy("$TARGET", "$SOURCE"))
+
+copy_files = [
+    ("tamaas/__init__.py", "#/python/tamaas.py"),
+    ("setup.py", "#/python/setup.py"),
+    ("tamaas/dumpers.py", "#/python/dumpers.py"),
+    ("tamaas/nonlinear_solvers.py", "#python/nonlinear_solvers.py")
+]
+
+for target, source in copy_files:
+    copy_env.Command(target, source, Copy("$TARGET", "$SOURCE"))
diff --git a/python/dumpers.py b/python/dumpers.py
new file mode 100644
index 0000000..4a5a55a
--- /dev/null
+++ b/python/dumpers.py
@@ -0,0 +1,187 @@
+# @author Lucas Frérot <lucas.frerot@epfl.ch>
+#
+# @section LICENSE
+#
+# Copyright (©)  2018 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/>.
+
+import numpy as np
+import scipy.interpolate as ip
+import tamaas as tm
+import os
+import sys
+
+try:
+    import uvw
+except ImportError:
+    print("uvw is not installed. You can install it with "
+          + "'pip install --user uvw'")
+    sys.exit(1)
+
+
+def makePeriodic(data):
+    data = np.append(data, np.expand_dims(data[:, 0, ...], axis=1), axis=1)
+    data = np.append(data, np.expand_dims(data[:, :, 0, ...], axis=2), axis=2)
+    return data
+
+
+def makeTractionPeriodic(data):
+    data = np.append(data, np.expand_dims(data[0, ...], axis=0), axis=0)
+    data = np.append(data, np.expand_dims(data[:, 0, ...], axis=1), axis=1)
+    return data
+
+
+class UVWDumper(tm.ModelDumper):
+    """Dumper for elasto-plastic calculations"""
+    def __init__(self, residual, filename, dump_numpys=False):
+        tm.ModelDumper.__init__(self)
+        self.residual = residual
+        self.count = 0
+        self.filename = filename
+        self.dump_numpys = dump_numpys
+
+        if not os.path.exists('paraview'):
+            os.mkdir('paraview')
+        if self.dump_numpys and not os.path.exists('numpys'):
+            os.mkdir('numpys')
+
+    def dump(self, model):
+        self.uvw_dump(model)
+        if self.dump_dumpys:
+            self.numpy_dump(model)
+        self.count += 1
+
+    def uvw_dump(self, model):
+        """Dump displacements, plastic deformations and stresses"""
+        discretization = model.getDiscretization().copy()
+        discretization[1] += 1
+        discretization[2] += 1
+
+        coordinates = [np.linspace(0, L, N)
+                       for L, N in zip(model.getSystemSize(),
+                                       discretization)]
+        # Correct order of coordinate dimensions
+        dimension_indices = [1, 2, 0]
+
+        # Creating rectilinear grid with correct order for components
+        grid = uvw.RectilinearGrid(
+            "paraview/{}_{:04d}.vtr".format(self.filename, self.count),
+            (coordinates[i] for i in dimension_indices))
+
+        # Point data arrays
+        stress = uvw.DataArray(makePeriodic(self.residual.getStress()),
+                               dimension_indices, 'stress')
+        plastic_strain = uvw.DataArray(
+            makePeriodic(self.residual.getPlasticStrain()),
+            dimension_indices, 'plastic_strain')
+        displacement = uvw.DataArray(makePeriodic(model.getDisplacement()),
+                                     dimension_indices, 'displacement')
+        residual = uvw.DataArray(makePeriodic(self.residual.getVector()),
+                                 dimension_indices, 'residual')
+
+        grid.addPointData(stress)
+        grid.addPointData(plastic_strain)
+        grid.addPointData(displacement)
+        grid.addPointData(residual)
+
+        grid.write()
+
+    def numpy_dump(self, model):
+        normal_traction = makeTractionPeriodic(model.getTraction()[:, :, 2])
+        normal_displacement = makeTractionPeriodic(
+            model.getDisplacement()[0, :, :, 2])
+        np.savez('numpys/bem_{:04d}.npz'.format(self.count),
+                 traction=normal_traction,
+                 displacement=normal_displacement)
+
+
+class AkantuDumper(tm.ModelDumper):
+    def __init__(self, ak_solver, filename, extra_dumps=False):
+        tm.ModelDumper.__init__(self)
+        self.solver = ak_solver
+        self.model = self.solver.model
+        self.filename = filename
+        self.model.setBaseName(filename)
+        self.extra_dumps = extra_dumps
+
+        if not os.path.exists('paraview'):
+            os.mkdir('paraview')
+        if self.extra_dumps and not os.path.exists('numpys'):
+            os.mkdir('numpys')
+
+        fields = ["displacement",
+                  "external_force",
+                  "internal_force",
+                  "inelastic_strain",
+                  "iso_hardening",
+                  "stress",
+                  "blocked_dofs"]
+
+        for field in fields:
+            self.model.addDumpField(field)
+
+        self.bot_face = np.where(self.model.getMesh().getNodes()[:, 2] == 0)
+        self.count = 0
+
+    def dump(self, model):
+        self.akantu_dump(model)
+        if self.extra_dumps:
+            self.uvw_dump(model)
+            self.numpy_dump(model)
+        self.count += 1
+
+    def uvw_dump(self, model):
+        self.discretization = np.array(model.getDiscretization()[1:], copy=True)
+        self.discretization += 1
+
+        self.coordinates = [np.linspace(0, L, N)
+                            for L, N in zip(model.getSystemSize()[1:],
+                                            self.discretization)]
+
+        # Correct order of coordinate dimensions
+        dimension_indices = range(2)
+
+        # Creating rectilinear grid with correct order for components
+        grid = uvw.RectilinearGrid(
+            "paraview/{}_traction_{:04d}.vtr".format(self.filename,
+                                                     self.count),
+            (self.coordinates[i] for i in dimension_indices))
+        traction = uvw.DataArray(makeTractionPeriodic(model.getTraction()),
+                                 dimension_indices, 'traction')
+        grid.addPointData(traction)
+        grid.write()
+
+    def akantu_dump(self, model):
+        self.model.dump()
+
+    def numpy_dump(self, model):
+        arrays = {}
+
+        arrays['traction'] = makeTractionPeriodic(model.getTraction()[:, :, 2])
+        displacements = self.model.getDisplacement()[self.bot_face, 2]
+
+        bot_points = self.model.getMesh().getNodes()[self.bot_face, :-1]
+        bot_points = bot_points.reshape(
+            (self.discretization[0] * self.discretization[1], 2))
+        x, y = np.meshgrid(*self.coordinates, indexing='ij')
+        disp = ip.griddata(bot_points,
+                           displacements.T, (x, y))
+        disp = disp.reshape(self.discretization)
+        arrays['displacement'] = disp - disp.mean()
+
+        np.savez('numpys/{}_{:04d}.npz'.format(self.filename,
+                                               self.count), **arrays)
diff --git a/python/nonlinear_solvers.py b/python/nonlinear_solvers.py
new file mode 100644
index 0000000..06efc23
--- /dev/null
+++ b/python/nonlinear_solvers.py
@@ -0,0 +1,316 @@
+# @author Lucas Frérot <lucas.frerot@epfl.ch>
+#
+# @section LICENSE
+#
+# Copyright (©)  2018 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/>.
+
+import tamaas as tm
+
+from scipy.optimize import root
+
+import numpy as np
+import scipy.interpolate as ip
+from scipy.sparse.linalg import LinearOperator, lgmres
+from scipy.linalg import norm
+from scipy.optimize import newton_krylov
+
+
+def makePeriodic(data):
+    data = np.append(data, np.expand_dims(data[0, ...], axis=0), axis=0)
+    data = np.append(data, np.expand_dims(data[:, 0, ...], axis=1), axis=1)
+    return data
+
+
+class NonSolver(tm.EPSolver):
+    def __init__(self, residual, model):
+        tm.EPSolver.__init__(self, residual)
+        self.model = model
+        self.elastic_strain_increment = self.getStrainIncrement()
+        self.last_increment = self.elastic_strain_increment.copy()
+
+    def solve(self):
+        self.elastic_strain_increment[...] = 0
+        self.getResidual().computeResidual(self.elastic_strain_increment)
+        self.elastic_strain_increment[...] = self.getResidual().getVector()
+        self.last_increment *= -1
+        self.last_increment += self.elastic_strain_increment
+        self.getResidual().computeResidualDisplacement(
+            self.elastic_strain_increment)
+
+
+class ScipySolver(tm.EPSolver):
+    def __init__(self, residual, model):
+        tm.EPSolver.__init__(self, residual)
+        self.model = model
+        self._x = self.getStrainIncrement()
+        self._residual = self.getResidual()
+        self.surface_correction = model.getDisplacement()[0, :, :, 2]
+        self.last_residual_disp = np.zeros_like(self.surface_correction)
+
+    def solve(self):
+        # 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.last_residual_disp[...] = 0
+        self._x[...] = 0
+
+
+class NewtonKrylovSolver(ScipySolver):
+    def __init__(self, residual, model):
+        ScipySolver.__init__(self, residual, model)
+        self.options = {'ftol': 1e-8, 'fatol': 1e-15}
+
+    def scipy_solve(self, compute_residual):
+        return newton_krylov(compute_residual, self._x,
+                             f_tol=1e-8, verbose=True)
+
+
+class DFSANESolver(ScipySolver):
+    def __init__(self, residual, model):
+        ScipySolver.__init__(self, residual, model)
+        self.options = {'ftol': 1e-10, 'fatol': 1e-10}
+
+    def scipy_solve(self, compute_residual):
+        solution = root(compute_residual, self._x, method='df-sane',
+                        options=self.options)
+        print("DF-SANE: {} ({} iterations)".format(solution.message,
+                                                   solution.nit))
+        return solution.x.copy()
+
+
+class NewtonRaphsonSolver(tm.EPSolver):
+    def __init__(self, residual, model):
+        tm.EPSolver.__init__(self, 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 FromTraction:
+    """Functor for Neumann boundary conditions"""
+    def __init__(self, interpolation):
+        self.interpolation = interpolation
+
+    def operator(self, quad_point, dual, coord, normals):
+        val = self.interpolation(*(coord[:-1]))
+        dual[2] = val[0]
+
+
+try:
+    import akantu as ak
+
+    class PoincareFEM_plast(tm.Residual):
+        def __init__(self, tm_model, ak_model):
+            tm.Residual.__init__(self, tm_model)
+            self.tm_model = tm_model
+            self.discretization = tm_model.getDiscretization().copy()
+            self.previous_pressure = np.zeros(self.discretization[1:])
+            self.residual_displacement = np.zeros(self.discretization[1:])
+            self.discretization[1] += 1
+            self.discretization[2] += 1
+            self.coordinates = [np.linspace(0, L, N)
+                                for L, N in zip(tm_model.getSystemSize(),
+                                                self.discretization)]
+
+            self.model = ak_model
+            self.bot_face = np.where(self.model.getMesh().getNodes()[:, 2] == 0)
+            self.ids = \
+                {'inelastic_strain':
+                 'previous_solid_mechanics_model:0:'
+                 + 'plastic_linear_isotropic_hardening:inelastic_strain',
+                 'iso_hardening':
+                 'previous_solid_mechanics_model:0:'
+                 + 'plastic_linear_isotropic_hardening:iso_hardening'}
+            self.savePreviousState()
+
+        def apply(self, field, out=None):
+            interpolation = ip.interp2d(self.coordinates[1],
+                                        self.coordinates[2],
+                                        makePeriodic(field))
+            functor = FromTraction(interpolation)
+            self.model.getExternalForce()[...] = 0
+            self.model.applyNeumannBC(functor, "bot")
+            self.model.solveStep()
+            displacements = self.model.getDisplacement()[self.bot_face, 2]
+
+            if out is not None:
+                # All these shenanigans because of node ordering
+                bot_points = self.model.getMesh().getNodes()[self.bot_face, :-1]
+                bot_points = bot_points.reshape(
+                    (self.discretization[1] * self.discretization[2], 2))
+                x, y = np.meshgrid(*self.coordinates[1:], indexing='ij')
+                disp = ip.griddata(bot_points, displacements.T, (x, y))
+                out[...] = disp.reshape(self.discretization[1:])[:-1, :-1]
+                out -= out.mean()
+
+        def savePreviousState(self):
+            self.inelastic_strain = \
+                self.model.getMaterial(0).getArrayReal(
+                    self.ids['inelastic_strain'], ak._tetrahedron_4).copy()
+            self.hardening = \
+                self.model.getMaterial(0).getArrayReal(
+                    self.ids['iso_hardening'], ak._tetrahedron_4).copy()
+            self.previous_pressure[...] = self.tm_model.getTraction()[:, :, 2]
+
+        def restorePreviousState(self):
+            inelastic_strain = \
+                self.model.getMaterial(0).getArrayReal(
+                    self.ids['inelastic_strain'], ak._tetrahedron_4)
+            hardening = \
+                self.model.getMaterial(0).getArrayReal(
+                    self.ids['iso_hardening'], ak._tetrahedron_4)
+            inelastic_strain[...] = self.inelastic_strain
+            hardening[...] = self.hardening
+
+        def updateState(self, increment):
+            """
+            Update the plastic deformation state
+            - Here we disregard the increment
+            """
+            self.restorePreviousState()
+            self.apply(self.tm_model.getTraction()[:, :, 2])
+            self.apply(np.zeros_like(self.residual_displacement),
+                       self.residual_displacement)
+            self.apply(self.tm_model.getTraction()[:, :, 2])
+            self.savePreviousState()
+
+        def getStrainIncrement(self):
+            return np.array()
+
+    class AkantuSolver(tm.EPSolver):
+        """Non-linear finite elements solver"""
+        def __init__(self, tm_model, material_file, mesh_name):
+            """Initialize Akantu SolidMechanicsModel"""
+            ak.parseInput(material_file)
+            self.tm_model = tm_model
+
+            self.name = mesh_name
+
+            # Loading mesh
+            self.mesh = ak.Mesh(3)
+            self.mesh.read(mesh_name)
+
+            # Creating Akantu model
+            self.model = ak.SolidMechanicsModel(self.mesh)
+            self.model.initFull(_analysis_method=ak._static)
+
+            # Poincare operator/dummy_residual
+            self.poincare = PoincareFEM_plast(tm_model, self.model)
+            tm.EPSolver.__init__(self, self.poincare)
+
+            # Setting boundary conditions
+            position = self.mesh.getNodes()
+            blocked_dofs = self.model.getBlockedDOFs()
+
+            epsilon = 1e-10
+            for i in range(self.mesh.getNbNodes()):
+                if np.abs(position[i, 2] - 1) < epsilon:  # Blocking top surface
+                    blocked_dofs[i, 2] = True
+                    if np.abs(position[i, 0]) < epsilon:
+                        blocked_dofs[i, 1] = True
+                        if np.abs(position[i, 1]) < epsilon:
+                            blocked_dofs[i, 0] = True
+                if np.abs(position[i, 0]) < epsilon \
+                   or np.abs(position[i, 0] - 1.) < epsilon:
+                    blocked_dofs[i, 0] = True
+                if np.abs(position[i, 1]) < epsilon \
+                   or np.abs(position[i, 1] - 1.) < epsilon:
+                    blocked_dofs[i, 1] = True
+
+            solver = self.model.getNonLinearSolver("static")
+            solver.set('max_iterations', 200)
+            solver.set('threshold', 1e-6)
+            solver.set('convergence_type', ak._scc_residual)
+
+            self.surface_correction = tm_model.getDisplacement()[0, :, :, 2]
+            self.previous_solution = np.zeros_like(self.surface_correction)
+
+        def solve(self):
+            """Apply and release pressure to get the residual displacements"""
+            # this is the pressure increment
+            pressure = self.tm_model.getTraction()[:, :, 2]
+            self.poincare.restorePreviousState()
+            self.poincare.apply(pressure + self.poincare.previous_pressure)
+            residual_disp = np.zeros_like(self.surface_correction)
+            self.poincare.apply(np.zeros_like(pressure), residual_disp)
+            self.surface_correction[...] = \
+                residual_disp - self.poincare.residual_displacement
+
+        def reset(self):
+            self.previous_solution[...] = 0
+
+except ImportError:
+    print("Akantu module not found: AkantuSolver unavailable")
diff --git a/python/setup.py b/python/setup.py
index 08b9e76..e3a6140 100644
--- a/python/setup.py
+++ b/python/setup.py
@@ -1,19 +1,19 @@
 import setuptools
 import sysconfig
 
 setuptools.setup(
     name="tamaas",
     version="1.0.0",
-    
+
     packages=setuptools.find_packages(),
     package_data={'tamaas': ['_tamaas.{}'.format(
         sysconfig.get_config_vars()['EXT_SUFFIX'])
     ]},
 
     author=['Lucas Frérot',
             'Guillaume Anciaux',
             'Valentine Rey',
             'Son Pham-Ba'],
     description='A fast rough contact library',
     url="https://c4science.ch/project/view/2036/",
 )