diff --git a/examples/.gitignore b/examples/.gitignore
index fb6ae0a..4a474f3 100644
--- a/examples/.gitignore
+++ b/examples/.gitignore
@@ -1,3 +1,4 @@
 paraview
 hdf5
 numpy*
+.gdb_history # yes I debug with examples ;-)
diff --git a/examples/plasticity.py b/examples/plasticity.py
index a86bc29..3896bee 100644
--- a/examples/plasticity.py
+++ b/examples/plasticity.py
@@ -1,79 +1,86 @@
 #!/usr/bin/env python3
 # @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/>.
 
 import numpy as np
 import tamaas as tm
 
-from tamaas.dumpers import UVWGroupDumper as Dumper
-from tamaas.nonlinear_solvers import DFSANESolver as Solver, ToleranceManager
+from tamaas.dumpers import H5Dumper as Dumper
+from tamaas.nonlinear_solvers import DFSANECXXSolver as Solver, ToleranceManager
 
 
+tm.initialize(1)
+
 # Definition of modeled domain
 model_type = tm.model_type.volume_2d
 discretization = [32, 51, 51]
 flat_domain = [1, 1]
 system_size = [0.5] + flat_domain
 
 # Creation of model
 model = tm.ModelFactory.createModel(model_type,
                                     system_size,
                                     discretization)
 
 model.E = 1.
 model.nu = 0.3
 
 # Setup for plasticity
 residual = tm.ModelFactory.createResidual(model,
                                           sigma_y=0.1 * model.E,
                                           hardening=0.01 * model.E)
 # Possibly change integration method
 # residual.setIntegrationMethod(tm.integration_method.cutoff, 1e-12)
 
 # Setup non-linear solver with variable tolerance
 epsolver = ToleranceManager(1e-5, 1e-9, 1/4)(Solver)(residual)
 
 # Setup for contact
 x = np.linspace(0, system_size[1], discretization[1],
                 endpoint=False, dtype=tm.dtype)
 y = np.linspace(0, system_size[2], discretization[2],
                 endpoint=False, dtype=tm.dtype)
 xx, yy = np.meshgrid(x, y, indexing='ij')
 
 R = 0.2
 surface = -((xx - flat_domain[0] / 2)**2 +
             (yy - flat_domain[1] / 2)**2) / (2 * R)
 
-csolver = tm.PolonskyKeerRey(model, surface, 1e-12,
+local_shape = tm.mpi.local_shape(surface.shape)
+offset = tm.mpi.local_offset(surface.shape)
+local_surface = surface[offset:offset+local_shape[0], :].copy()
+
+csolver = tm.PolonskyKeerRey(model, local_surface, 1e-12,
                              tm.PolonskyKeerRey.pressure,
                              tm.PolonskyKeerRey.pressure)
 
 # EPIC setup
 epic = tm.EPICSolver(csolver, epsolver, 1e-7)
 
 # Dumper
 dumper_helper = Dumper('hertz', 'displacement', 'stress', 'plastic_strain')
 model.addDumper(dumper_helper)
 
 loads = np.linspace(0.001, 0.005, 3)
 
 for i, load in enumerate(loads):
     epic.acceleratedSolve(load)
     model.dump()
-    print("---> Solved load step {}/{}".format(i+1, len(loads)))
+    tm.Logger().get(tm.LogLevel.info) \
+        << "---> Solved load step {}/{}".format(i+1, len(loads))