diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp
index acbc22b..129ab55 100644
--- a/src/core/statistics.cpp
+++ b/src/core/statistics.cpp
@@ -1,214 +1,214 @@
 /**
  *  @file
  *  LICENSE
  *
  *  Copyright (©) 2016-2021 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 "statistics.hh"
 #include "fft_engine.hh"
 #include "loop.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 template <UInt dim>
 Real Statistics<dim>::computeRMSHeights(Grid<Real, dim>& surface) {
   return std::sqrt(surface.var());
 }
 
 template <UInt dim>
 Real Statistics<dim>::computeSpectralRMSSlope(Grid<Real, dim>& surface) {
   const auto h_size =
       GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(h_size);
   wavevectors *= 2 * M_PI;  // need q for slopes
   const auto psd = computePowerSpectrum(surface);
 
   const Real rms_slope_mean = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(VectorProxy<const Real, dim> q, const Complex& psd_val) {
         // Checking if we're in the zone that does not have hermitian symmetry
         if (std::abs(q.back()) < 1e-15)
           return q.l2squared() * psd_val.real();
         else
           return 2 * q.l2squared() * psd_val.real();
       },
       range<VectorProxy<const Real, dim>>(wavevectors), psd);
 
   return std::sqrt(rms_slope_mean);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 GridHermitian<Real, dim>
 Statistics<dim>::computePowerSpectrum(Grid<Real, dim>& surface) {
   const auto h_size =
       GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
   GridHermitian<Real, dim> psd(h_size, surface.getNbComponents());
 
   FFTEngine::makeEngine()->forward(surface, psd);
 
   Real factor = 1. / surface.getGlobalNbPoints();
 
   // Squaring the fourier transform of surface and normalizing
   Loop::loop(
       [factor] CUDA_LAMBDA(Complex & c) {
         c *= factor;
         c *= conj(c);
       },
       psd);
 
   return psd;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 Grid<Real, dim>
 Statistics<dim>::computeAutocorrelation(Grid<Real, dim>& surface) {
   Grid<Real, dim> acf(surface.sizes(), surface.getNbComponents());
   auto psd = computePowerSpectrum(surface);
 
   FFTEngine::makeEngine()->backward(acf, psd);
   acf *= acf.getGlobalNbPoints();
   return acf;
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 Real Statistics<dim>::contact(const Grid<Real, dim>& tractions,
                               UInt perimeter) {
   Real points = 0;
   UInt nc = tractions.getNbComponents();
 
   switch (nc) {
   case 1:
     points = Loop::reduce(
         [] CUDA_LAMBDA(const Real& t) -> Real { return t > 0; }, tractions);
     break;
   case 2:
     points = Loop::reduce(
         [] CUDA_LAMBDA(VectorProxy<const Real, 2> t) -> Real {
           return t.back() > 0;
         },
         range<VectorProxy<const Real, 2>>(tractions));
     break;
   case 3:
     points = Loop::reduce(
         [] CUDA_LAMBDA(VectorProxy<const Real, 3> t) -> Real {
           return t.back() > 0;
         },
         range<VectorProxy<const Real, 3>>(tractions));
     break;
   default:
     TAMAAS_EXCEPTION("Invalid number of components in traction");
   }
 
-  auto area = points / tractions.getNbPoints();
+  auto area = points / tractions.getGlobalNbPoints();
 
   if (dim == 1)
     perimeter = 0;
 
   // Correction from Yastrebov et al. (Trib. Intl., 2017)
   // 10.1016/j.triboint.2017.04.023
   return area -
-         (M_PI - 1 + std::log(2)) / (24. * tractions.getNbPoints()) * perimeter;
+         (M_PI - 1 + std::log(2)) / (24. * tractions.getGlobalNbPoints()) * perimeter;
 }
 
 /* -------------------------------------------------------------------------- */
 namespace {
 template <UInt dim>
 class moment_helper {
 public:
   moment_helper(const std::array<UInt, dim>& exp) : exponent(exp) {}
 
   CUDA_LAMBDA Complex operator()(VectorProxy<Real, dim> q,
                                  const Complex& phi) const {
     Real mul = 1;
     for (UInt i = 0; i < dim; ++i)
       mul *= std::pow(q(i), exponent[i]);
 
     // Do not duplicate everything from hermitian symmetry
     if (std::abs(q.back()) < 1e-15)
       return mul * phi;
     else
       return 2 * mul * phi;
   }
 
 private:
   std::array<UInt, dim> exponent;
 };
 }  // namespace
 
 template <>
 std::vector<Real> Statistics<1>::computeMoments(Grid<Real, 1>& surface) {
   constexpr UInt dim = 1;
   std::vector<Real> moments(3);
   const auto psd = computePowerSpectrum(surface);
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(psd.sizes());
   // we don't multiply by 2 pi because moments are computed with k
   moments[0] = Loop::reduce<operation::plus>(moment_helper<dim>{{{0}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   moments[1] = Loop::reduce<operation::plus>(moment_helper<dim>{{{2}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   moments[2] = Loop::reduce<operation::plus>(moment_helper<dim>{{{4}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   return moments;
 }
 
 template <>
 std::vector<Real> Statistics<2>::computeMoments(Grid<Real, 2>& surface) {
   constexpr UInt dim = 2;
   std::vector<Real> moments(3);
   const auto psd = computePowerSpectrum(surface);
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(psd.sizes());
   // we don't multiply by 2 pi because moments are computed with k
   moments[0] = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 0}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   auto m02 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 2}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m20 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 0}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   moments[1] = 0.5 * (m02 + m20);
   auto m22 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 2}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m40 = Loop::reduce<operation::plus>(moment_helper<dim>{{{4, 0}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m04 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 4}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   moments[2] = (3 * m22 + m40 + m04) / 3.;
   return moments;
 }
 
 template struct Statistics<1>;
 template struct Statistics<2>;
 
 }  // namespace tamaas
diff --git a/tests/mpi_routines.py b/tests/mpi_routines.py
index 33cb5ec..0817740 100644
--- a/tests/mpi_routines.py
+++ b/tests/mpi_routines.py
@@ -1,281 +1,284 @@
 # -*- coding: utf-8 -*-
 # @file
 # LICENSE
 #
 # Copyright (©) 2016-2021 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 print_function
 import shutil
 
 import tamaas as tm
 import numpy as np
 
 from mpi4py import MPI
 from numpy.linalg import norm
 
 from conftest import HertzFixture
 
 
 def make_surface(N):
     spectrum = tm.Isopowerlaw2D()
     spectrum.q0 = 2
     spectrum.q1 = 2
     spectrum.q2 = 16
     spectrum.hurst = 0.8
 
     generator = tm.SurfaceGeneratorRandomPhase2D(N)
     generator.spectrum = spectrum
     generator.random_seed = 0
 
     return generator.buildSurface()
 
 
 def mpi_surface_generator():
     N = [128, 128]
     tm.set_log_level(tm.LogLevel.debug)
 
     seq_surface = None
     comm = MPI.COMM_WORLD
     print('[{}] {}'.format(comm.rank, comm.size))
 
     with tm.mpi.sequential():
         if comm.rank == 0:
             seq_surface = make_surface(N)
 
     surface = make_surface(N)
     print('[{}] {}'.format(comm.rank, surface.shape))
 
     recv = comm.gather(surface, root=0)
 
     if comm.rank == 0:
         gsurface = np.concatenate(recv)
 
         if False:
             import matplotlib.pyplot as plt
             plt.imshow(seq_surface)
             plt.colorbar()
             plt.figure()
             plt.imshow(gsurface)
             plt.colorbar()
             plt.show()
 
         # I think the only reason this assert works is because the
         # spectrum is compact in the Fourier domain -> surface is regular
         assert np.all(seq_surface == gsurface)
 
 
 def mpi_model_creation():
     N = [20, 50, 50]
     S = [1., 1., 1.]
     comm = MPI.COMM_WORLD
 
     def get_discretizations(model):
         return model.shape, model.boundary_shape
 
     model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S[1:], N[1:])
     n, bn = get_discretizations(model)
     n[0] = comm.allreduce(n[0])
     bn[0] = comm.allreduce(bn[0])
     assert n == N[1:] and bn == N[1:]
 
     model = tm.ModelFactory.createModel(tm.model_type.volume_2d, S, N)
     n, bn = get_discretizations(model)
     n[1] = comm.allreduce(n[1])
     bn[0] = comm.allreduce(bn[0])
     assert n == N and bn == N[1:]
 
 
 def mpi_polonsky_keer():
     N = [1024, 1024]
     S = [1., 1.]
     load = 0.00001
     comm = MPI.COMM_WORLD
 
     hertz = HertzFixture(N[0], load)
     surface = hertz.surface
 
     tm.set_log_level(tm.LogLevel.debug)
 
     def solve():
         model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N)
         model.E, model.nu = 1, 0
 
         local_n, local_offset = tm.mpi.local_shape(N), tm.mpi.local_offset(N)
         local_surface = np.copy(
             surface[local_offset:local_offset+local_n[0], :])
 
         solver = tm.PolonskyKeerRey(model, local_surface, 1e-15)
 
         print('Created solver')
         solver.solve(load)
         return np.copy(model['traction']), np.copy(model['displacement'])
 
     tractions, displacements = solve()
 
     precv = comm.gather(tractions, root=0)
     drecv = comm.gather(displacements, root=0)
 
     if comm.rank == 0:
         tractions = np.concatenate(precv)
         displacements = np.concatenate(drecv)
 
         perror = norm(tractions - hertz.pressure) / norm(hertz.pressure)
         derror = norm(displacements - hertz.displacement)\
             / norm(hertz.displacement)
 
         if False:
             print(perror, derror)
             import matplotlib.pyplot as plt
             plt.imshow(tractions - hertz.pressure)
             plt.colorbar()
             plt.show()
 
         assert perror < 5e-3
         assert derror < 8e-3
 
 
 def mpi_polonsky_keer_compare():
     N = [273, 273]
     S = [1., 1.]
     E, nu = 0.2, 0.3
     load = 0.1
     comm = MPI.COMM_WORLD
 
     seq_tractions = None
     rms = 0
 
     with tm.mpi.sequential():
         if comm.rank == 0:
             model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N)
             model.E, model.nu = E, nu
             surface = make_surface(N)
             rms = tm.Statistics2D.computeSpectralRMSSlope(surface)
             surface /= rms
             solver = tm.PolonskyKeerRey(model, surface, 1e-15)
             solver.solve(load)
-            seq_tractions = np.copy(model['traction'])
+            seq_tractions = np.copy(model.traction)
+            seq_area = tm.Statistics2D.contact(model.traction)
 
     rms = comm.bcast(rms, root=0)
 
     model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N)
     model.E, model.nu = E, nu
     surface = make_surface(N) / rms
     solver = tm.PolonskyKeerRey(model, surface, 1e-15)
     solver.solve(load)
     tractions = model['traction']
+    area = tm.Statistics2D.contact(tractions)
 
     recv = comm.gather(tractions, root=0)
 
     if comm.rank == 0:
         tractions = np.concatenate(recv)
         error = np.linalg.norm(seq_tractions - tractions) / seq_tractions.size
 
         if False:
             print(error)
             import matplotlib.pyplot as plt
             plt.imshow(seq_tractions - tractions)
             plt.colorbar()
             plt.show()
 
         assert error < 1e-13
+        assert np.abs(seq_area - area) / seq_area < 1e-13
 
 
 def mpi_h5dump():
     try:
         from tamaas.dumpers import H5Dumper as Dumper
         import h5py
     except ImportError:
         return
 
     if not h5py.get_config().mpi:
         print('Did not test H5Dumper with MPI')
         return
 
     model_type = tm.model_type.volume_2d
     discretization = [10, 2, 10]
     flat_domain = [1, 1]
     system_size = [0.5] + flat_domain
 
     # Creation of model
     model = tm.ModelFactory.createModel(model_type,
                                         system_size,
                                         discretization)
 
     model['displacement'][:] = MPI.COMM_WORLD.rank
     model['traction'][:] = MPI.COMM_WORLD.rank
 
     dumper_helper = Dumper('test_mpi_dump', 'displacement', 'traction')
     dumper_helper << model
 
     with h5py.File('hdf5/test_mpi_dump_0000.h5', 'r', driver='mpio',
                    comm=MPI.COMM_WORLD) as handle:
         assert np.all(handle['displacement'][:, 0, :] == 0)
         assert np.all(handle['displacement'][:, 1, :] == 1)
         assert np.all(handle['traction'][0, :] == 0)
         assert np.all(handle['traction'][1, :] == 1)
 
     rmodel = dumper_helper.read('hdf5/test_mpi_dump_0000.h5')
 
     assert np.all(rmodel.displacement == MPI.COMM_WORLD.rank)
     assert np.all(rmodel.traction == MPI.COMM_WORLD.rank)
 
     if tm.mpi.rank() == 0:
         shutil.rmtree('hdf5')
 
 
 def mpi_plastic_solve():
     from conftest import UniformPlasticity
     from tamaas.nonlinear_solvers import DFSANECXXSolver as Solver
 
     patch_isotropic_plasticity = UniformPlasticity(tm.model_type.volume_2d,
                                                    [1.] * 3,
                                                    [4] * 3)
 
     model = patch_isotropic_plasticity.model
     residual = patch_isotropic_plasticity.residual
 
     applied_pressure = 0.1
 
     solver = Solver(residual)
     solver.tolerance = 1e-15
     pressure = model['traction'][..., 2]
     pressure[:] = applied_pressure
 
     solver.solve()
     solver.updateState()
 
     solution, normal = patch_isotropic_plasticity.solution(applied_pressure)
 
     for key in solution:
         error = norm(model[key] - solution[key]) / normal[key]
         assert error < 2e-15
 
 
 def mpi_flood_fill():
     contact = np.zeros([10, 1], dtype=np.bool)
     contact[:1, :] = True
     contact[-1:, :] = True
 
     clusters = tm.FloodFill.getClusters(contact, False)
 
     if tm.mpi.rank() == 0:
         assert len(clusters) == 3, "Wrong number of clusters"
         assert [c.getArea() for c in clusters] == [1, 2, 1], "Wrong areas"
 
 
 if __name__ == '__main__':
     mpi_polonsky_keer_compare()