diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp
index 8f362f5..edc1cab 100644
--- a/python/wrap/solvers.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,119 +1,120 @@
 /**
  * @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 "wrap.hh"
 #include "contact_solver.hh"
 #include "polonsky_keer_rey.hh"
 #include "kato.hh"
 #include "beck_teboulle.hh"
 #include "condat.hh"
 #include "polonsky_keer_tan.hh"
 #include "epic.hh"
 #include "ep_solver.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 namespace wrap {
 
 /* -------------------------------------------------------------------------- */
 using namespace py::literals;
 
 class PyEPSolver : public EPSolver {
 public:
   using EPSolver::EPSolver;
   void solve() override {
     PYBIND11_OVERLOAD_PURE(void, EPSolver, solve);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 void wrapSolvers(py::module& mod) {
   py::class_<ContactSolver>(mod, "ContactSolver")
       // .def(py::init<Model&, const GridBase<Real>&, Real>(), "model"_a,
       //      "surface"_a, "tolerance"_a)
       .def("setMaxIterations", &ContactSolver::setMaxIterations, "max_iter"_a)
       .def("setDumpFrequency", &ContactSolver::setDumpFrequency, "dump_freq"_a)
       .def("addFunctionalTerm", &ContactSolver::addFunctionalTerm);
 
   py::class_<PolonskyKeerRey, ContactSolver> pkr(mod, "PolonskyKeerRey");
   pkr.def(py::init<Model&, const GridBase<Real>&, Real, PolonskyKeerRey::type,
                    PolonskyKeerRey::type>(),
           "model"_a, "surface"_a, "tolerance"_a, "primal_type"_a,
           "constraint_type"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
       .def("solve", &PolonskyKeerRey::solve, "target"_a)
       .def("computeError", &PolonskyKeerRey::computeError);
 
   py::enum_<PolonskyKeerRey::type>(pkr, "type")
       .value("gap", PolonskyKeerRey::gap)
       .value("pressure", PolonskyKeerRey::pressure)
       .export_values();
 
   py::class_<Kato, ContactSolver> kato(mod, "Kato");
   kato.def(py::init<Model&, const GridBase<Real>&, Real, Real>(),
           "model"_a, "surface"_a, "tolerance"_a, "mu"_a,
           py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
       .def("solve", &Kato::solve, "p0"_a, "proj_iter"_a=50)
       .def("solveRelaxed", &Kato::solveRelaxed, "g0"_a)
       .def("solveRegularized", &Kato::solveRegularized, "p0"_a, "r"_a=0.01)
       .def("computeCost", &Kato::computeCost, "use_tresca"_a=false);
 
   py::class_<BeckTeboulle, ContactSolver> bt(mod, "BeckTeboulle");
   bt.def(py::init<Model&, const GridBase<Real>&, Real, Real>(),
           "model"_a, "surface"_a, "tolerance"_a, "mu"_a,
           py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
       .def("solve", &BeckTeboulle::solve, "p0"_a)
       .def("computeCost", &BeckTeboulle::computeCost);
 
   py::class_<Condat, ContactSolver> cd(mod, "Condat");
   cd.def(py::init<Model&, const GridBase<Real>&, Real, Real>(),
           "model"_a, "surface"_a, "tolerance"_a, "mu"_a,
           py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
       .def("solve", &Condat::solve, "p0"_a, "grad_step"_a=0.9)
       .def("computeCost", &Condat::computeCost);
 
   py::class_<PolonskyKeerTan, ContactSolver> pkt(mod, "PolonskyKeerTan");
   pkt.def(py::init<Model&, const GridBase<Real>&, Real, Real>(),
           "model"_a, "surface"_a, "tolerance"_a, "mu"_a,
           py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
       .def("solve", &PolonskyKeerTan::solve, "p0"_a)
       .def("solveTresca", &PolonskyKeerTan::solveTresca, "p0"_a)
       .def("computeCost", &PolonskyKeerTan::computeCost, "use_tresca"_a=false);
   py::class_<EPSolver, PyEPSolver>(mod, "EPSolver")
       .def(py::init<Residual&>())
       .def("solve", &EPSolver::solve)
       .def("getStrainIncrement", &EPSolver::getStrainIncrement)
-      .def("getResidual", &EPSolver::getResidual);
+      .def("getResidual", &EPSolver::getResidual)
+      .def("updateState", &EPSolver::updateState);
 
   py::class_<EPICSolver>(mod, "EPICSolver")
       .def(py::init<ContactSolver&, EPSolver&, Real>(), "contact_solver"_a,
            "elasto_plastic_solver"_a, "tolerance"_a = 1e-10)
       .def("solve", &EPICSolver::solve, "load"_a);
 }
 
 }  // namespace wrap
 
 __END_TAMAAS__
diff --git a/src/SConscript b/src/SConscript
index 3ff0802..b005dcf 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,129 +1,131 @@
 import os
 Import('main_env')
 
 
 def prepend(path, list):
     return [os.path.join(path, x) for x in list]
 
 
-env = main_env
+env = main_env.Clone()
 
 # 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
 volume_potential.cpp
 kelvin.cpp
 mindlin.cpp
 boussinesq.cpp
 
 elasto_plastic/isotropic_hardening.cpp
 elasto_plastic/elasto_plastic_model.cpp
 elasto_plastic/residual.cpp
 """.split()
 model_list = prepend('model', model_list)
 
 # Solvers
 solvers_list = """
 contact_solver.cpp
 polonsky_keer_rey.cpp
 kato.cpp
 beck_teboulle.cpp
 condat.cpp
 polonsky_keer_tan.cpp
 ep_solver.cpp
 epic.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
 
 # For some reason collapse OMP loops don't compile on intel
 if env['CXX'] != 'icpc':
     rough_contact_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)
 
+# Adding extra warnings for Tamaas base lib
+env.AppendUnique(CXXFLAGS=['-Wextra'])
 libTamaas = env.SharedLibrary('Tamaas', rough_contact_list)
 Export('libTamaas')
diff --git a/src/bem/bem_functional.cpp b/src/bem/bem_functional.cpp
index 4477583..a71e9c3 100644
--- a/src/bem/bem_functional.cpp
+++ b/src/bem/bem_functional.cpp
@@ -1,44 +1,44 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 "bem_functional.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 Functional::Functional(BemFFTBase& bem) : parameters(), bem(bem), gradF(nullptr) {}
 
 Functional::~Functional() {}
 
 void Functional::setGradFPointer(Surface<Real>* p) { this->gradF = p; }
 
-void Functional::copyParameters(const std::map<std::string, Real>& params) {
+void Functional::copyParameters(const std::map<std::string, Real>& /*params*/) {
   // auto it = params.begin();
   // auto end = params.end();
   // for (; it != end ; ++it) this->parameters.insert(*it);
 }
 
 __END_TAMAAS__
diff --git a/src/bem/bem_gigi.cpp b/src/bem/bem_gigi.cpp
index 49f8e89..c55cd03 100644
--- a/src/bem/bem_gigi.cpp
+++ b/src/bem/bem_gigi.cpp
@@ -1,280 +1,280 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 "bem_gigi.hh"
 #include "surface.hh"
 #include <fstream>
 #include <iomanip>
 #include <iostream>
 #include <sstream>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 Real BemGigi::computeEquilibrium(Real epsilon, Real mean_displacement) {
 
   this->computeSpectralInfluenceOverDisplacement();
 
   Real Rold = 1.;
   Real f = 1e300;
   this->search_direction = 0.;
   this->true_displacements = 0;
   this->true_displacements = mean_displacement;
   this->computeGaps();
   this->optimizeToMeanDisplacement(mean_displacement);
   this->computeGaps();
 
   convergence_iterations.clear();
   nb_iterations = 0;
   max_iterations = 50;
   Real R = 0;
 
   while (f > epsilon && nb_iterations++ < max_iterations) {
 
     this->computeGaps();
     this->functional->computeGradFU();
     this->updateSearchDirection(R / Rold);
 
     Rold = R;
     Real alpha = this->computeOptimalStep();
     std::cout << "alpha vaut " << alpha << std::endl;
 
     this->old_displacements = this->true_displacements;
     this->updateDisplacements(alpha);
     this->computeGaps();
 
     this->optimizeToMeanDisplacement(mean_displacement);  // espace admissible
     this->computeGaps();
     this->computePressures(mean_displacement);
     f = computeStoppingCriterion();
   }
   this->computePressures(mean_displacement);
 
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeStoppingCriterion() {
 
   Real rho = this->functional->getParameter("rho");
   Real surface_energy = this->functional->getParameter("surface_energy");
 
   Real crit = 0.;
   // Real disp_norm = 0.;
 
   UInt n = surface.size();
   UInt size = n * n;
 
   //#pragma omp parallel for reduction(+:crit, disp_norm)
 #pragma omp parallel for reduction(+ : crit)
   for (UInt i = 0; i < size; ++i) {
     crit +=
         std::abs(this->gap(i) * (this->surface_tractions(i) +
                                  surface_energy / rho * exp(-(gap(i)) / rho)));
   }
 
   return crit;
 }
 
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::optimizeToMeanDisplacement(Real imposed_mean) {
   Real target_value =
       imposed_mean - SurfaceStatistics::computeAverage(surface, 0);
   UInt n = surface.size();
   UInt size = n * n;
 
   // Initial guesses for upper and lower bound
   Real step_min = -10;
   Real step_max = 10;
 
   // Gaps for upper and lower bound
   Real gap_min = positiveGapAverage(step_min);
   Real gap_max = positiveGapAverage(step_max);
 
   UInt max_expansion = 8;
 
   // Expand bounds if necessary
   for (UInt i = 0; gap_max < target_value && i < max_expansion;
        i++, step_max *= 10)
     gap_max = positiveGapAverage(step_max);
 
   for (UInt i = 0; gap_min > target_value && i < max_expansion;
        i++, step_min *= 10)
     gap_min = positiveGapAverage(step_min);
 
   Real g = 0.;
   Real epsilon = 1e-12;
   Real step = 0;
 
   while (fabs(g - target_value) > epsilon) {
     step = (step_min + step_max) / 2.;
     g = positiveGapAverage(step);
     if (g > target_value)
       step_max = step;
     else if (g < target_value)
       step_min = step;
     else {
       step_max = step;
       step_min = step;
     }
   }
 
   step = (step_min + step_max) / 2.;
 
 #pragma omp parallel for
   for (UInt i = 0; i < size; i++) {
     gap(i) += step;
     if (gap(i) < 0)
       gap(i) = 0;
     true_displacements(i) = gap(i) + surface(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::positiveGapAverage(Real shift) {
   UInt n = surface.size();
 
   Real res = 0;
 #pragma omp parallel for reduction(+ : res)
   for (UInt i = 0; i < n * n; i++) {
     Real shifted_gap = gap(i) + shift;
     res += shifted_gap * (shifted_gap > 0);
   }
 
   return res / (n * n);
 }
 
 /* -------------------------------------------------------------------------- */
 
-void BemGigi::updateSearchDirection(Real factor) {
+void BemGigi::updateSearchDirection(Real /*factor*/) {
   UInt n = surface.size();
   UInt size = n * n;
   const Surface<Real>& gradF = this->functional->getGradF();
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
 
     this->search_direction(i) = gradF(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeOptimalStep() {
   this->applyInverseInfluenceFunctions(search_direction, projected_direction);
 
   UInt n = surface.size();
   UInt size = n * n;
 
   const Surface<Real>& gradF = this->functional->getGradF();
   Real numerator = 0., denominator = 0.;
 
 #pragma omp parallel for reduction(+ : numerator, denominator)
   for (UInt i = 0; i < size; ++i) {
     numerator += gradF(i) * search_direction(i);
     denominator += projected_direction(i) * search_direction(i);
   }
 
   Real alpha = numerator / denominator;
   return alpha;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeTau() { return computeOptimalStep(); }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::updateDisplacements(Real alpha) {
   UInt n = surface.size();
   UInt size = n * n;
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     this->true_displacements(i) -= alpha * this->search_direction(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::emptyOverlap() {
   UInt n = surface.size();
   UInt size = n * n;
 
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     if (gap(i) < 0)
       this->true_displacements(i) = this->surface(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::enforceMeanDisplacement(Real mean_displacement) {
   UInt n = surface.size();
   UInt size = n * n;
   Real moyenne_surface = SurfaceStatistics::computeAverage(this->surface, 0);
   Real average = SurfaceStatistics::computeAverage(this->true_displacements, 0);
   Real factor =
       (mean_displacement - moyenne_surface) / (average - moyenne_surface);
 
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     this->true_displacements(i) =
         factor * (this->true_displacements(i) - this->surface(i)) +
         this->surface(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
-void BemGigi::computePressures(Real mean_displacement) {
+void BemGigi::computePressures(Real /*mean_displacement*/) {
   this->computeTractionsFromDisplacements();
 
   UInt n = surface.size();
   UInt size = n * n;
   Real rho = this->functional->getParameter("rho");
   Real surface_energy = this->functional->getParameter("surface_energy");
 
   Real mini = 3000;
   for (UInt i = 0; i < size; ++i) {
     if (gap(i) < rho) {
       if (this->surface_tractions(i) + surface_energy / rho < mini)
         mini = this->surface_tractions(i) + surface_energy / rho;
     } else {
       if (this->surface_tractions(i) < mini)
         mini = this->surface_tractions(i);
     }
   }
 
   this->surface_tractions -= mini;
 }
 
 __END_TAMAAS__
diff --git a/src/bem/bem_gigipol.cpp b/src/bem/bem_gigipol.cpp
index 551a5b6..00356b3 100644
--- a/src/bem/bem_gigipol.cpp
+++ b/src/bem/bem_gigipol.cpp
@@ -1,755 +1,755 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 "bem_gigipol.hh"
 #include "surface.hh"
 #include <fstream>
 #include <iomanip>
 #include <iostream>
 #include <sstream>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 void BemGigipol::computeTractionsFromDisplacements() {
   this->applyInverseInfluenceFunctions(this->true_displacements,
                                        this->surface_tractions);
 }
 
 Real BemGigipol::computeEquilibrium(Real epsilon, Real mean_displacement) {
 
   this->computeSpectralInfluenceOverDisplacement();
   this->surface_t = 0.;
   this->surface_r = 0.;
   this->search_direction = 0.;
   this->pold = 0.;
   this->surface_displacements = 0.;
   Real delta = 0.;
   Real Gold = 1.;
   Real f = 1e300;
   Real fPrevious = 1e300;
   UInt n = surface.size();
   UInt size = n * n;
 
   Real current_disp = SurfaceStatistics::computeAverage(true_displacements, 0);
   std::cout << "current disp " << current_disp << std::endl;
   if (current_disp <= 0.) {
     true_displacements = mean_displacement;
     std::cout << "je re initialise " << std::endl;
   }
 
   this->computeGaps();
   this->optimizeToMeanDisplacement(mean_displacement);
 
   std::ofstream file("output.txt");
   this->computeGaps();
 
   convergence_iterations.clear();
   nb_iterations = 0;
   max_iterations = 10000;
 
   while (f > epsilon && nb_iterations < max_iterations) {
     fPrevious = f;
 
     this->functional->computeGradFU();
     this->search_direction = this->functional->getGradF();
     Real gbar = this->computeMeanPressuresInNonContact();
 
     this->search_direction -= gbar;
     Real G = this->computeG();
     this->updateT(G, Gold, delta);
     Real tau = this->computeTau();
 
     this->old_displacements = this->true_displacements;
     Gold = G;
 #pragma omp parallel for
     for (unsigned int i = 0; i < size; ++i) {
 
       this->true_displacements(i) -= tau * this->surface_t(i);
     }
 
     // Projection on admissible space
     this->computeGaps();
 
 #pragma omp parallel for
     for (unsigned int i = 0; i < size; ++i) {
       if (this->gap(i) < 0) {
         this->true_displacements(i) = this->surface(i);
       }
     }
     this->computeGaps();
     delta = this->updateDisplacements(tau, mean_displacement);
 
     // Projection on admissible space
     this->computeGaps();
     this->enforceMeanDisplacement(mean_displacement);
 
     this->computeGaps();
     this->computePressures();
     f = this->computeStoppingCriterion();
 
     if (nb_iterations % dump_freq == 0) {
 
       std::cout << "G vaut " << G << std::endl;
       std::cout << "f vaut " << f << std::endl;
       std::cout << std::scientific << std::setprecision(10) << nb_iterations
                 << " " << f << " " << f - fPrevious << " " << 'G' << " " << G
                 << " " << std::endl;
 
       UInt n = surface.size();
       UInt size = n * n;
       Real crit = 0;
       Real disp_norm = 0;
       for (UInt i = 0; i < size; ++i) {
 
         crit += this->search_direction(i) * this->search_direction(i);
         disp_norm += (true_displacements(i) * true_displacements(i));
       }
 
       file << std::scientific << std::setprecision(10) << nb_iterations << " "
            << crit / disp_norm << " " << f << std::endl;
     }
 
     convergence_iterations.push_back(f);
     ++nb_iterations;
   }
 
   this->computePressures();
 
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 Real BemGigipol::computeEquilibrium2(Real epsilon, Real mean_pressure) {
 
   this->computeSpectralInfluenceOverDisplacement();
   this->surface_t = 0.;
   this->surface_r = 0.;
   this->search_direction = 0.;
   this->pold = 0.;
   this->surface_displacements = 0.;
   Real delta = 0.;
   Real Gold = 1.;
   Real f = 1e300;
   Real fPrevious = 1e300;
 
   UInt n = surface.size();
   UInt size = n * n;
 
   Real current_disp = SurfaceStatistics::computeAverage(surface, 0.);
   true_displacements = current_disp;
 
   this->computeGaps();
 
 #pragma omp parallel for
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) < 0) {
       this->true_displacements(i) = this->surface(i);
     }
   }
   this->computeGaps();
 
   // this->optimizeToMeanDisplacement(mean_displacement);
 
   std::ofstream file("output.txt");
 
   convergence_iterations.clear();
   nb_iterations = 0;
   max_iterations = 200000;
 
   while (f > epsilon && nb_iterations < max_iterations) {
     fPrevious = f;
 
     this->functional->computeGradFU();
     this->search_direction = this->functional->getGradF();
     Real gbar = this->computeMeanPressuresInNonContact();
 
     this->search_direction += (2 * mean_pressure + gbar);
 
     Real G = this->computeG();
 
     this->updateT(G, Gold, delta);
     Real tau = this->computeTau(mean_pressure);
     this->old_displacements = this->true_displacements;
     Gold = G;
 #pragma omp parallel for
     for (unsigned int i = 0; i < size; ++i) {
 
       this->true_displacements(i) -= tau * this->surface_t(i);
     }
 
     // Projection on admissible space
     this->computeGaps();
 
 #pragma omp parallel for
     for (unsigned int i = 0; i < size; ++i) {
       if (this->gap(i) < 0) {
         this->true_displacements(i) = this->surface(i);
       }
     }
     this->computeGaps();
     delta = this->updateDisplacements(tau, 0.);
 
     // Projection on admissible space
     this->computeGaps();
 
     this->computePressures();
     f = this->computeStoppingCriterion();
 
     this->computeTractionsFromDisplacements();
     this->functional->computeGradFU();
     const Surface<Real>& gradF = this->functional->getGradF();
     Real min = SurfaceStatistics::computeMinimum(gradF);
     if (f < epsilon && epsilon < std::abs(min + mean_pressure))
       f = 3.;
 
     if (nb_iterations % dump_freq == 0) {
 
       std::cout << "G vaut " << G << std::endl;
       std::cout << "f vaut " << f << std::endl;
       std::cout << std::scientific << std::setprecision(10) << nb_iterations
                 << " " << f << " " << f - fPrevious << " " << 'G' << " " << G
                 << " " << std::endl;
 
       UInt n = surface.size();
       UInt size = n * n;
       Real crit = 0;
       Real disp_norm = 0;
       for (UInt i = 0; i < size; ++i) {
 
         crit += this->search_direction(i) * this->search_direction(i);
         disp_norm += (true_displacements(i) * true_displacements(i));
       }
 
       file << std::scientific << std::setprecision(10) << nb_iterations << " "
            << crit / disp_norm << " " << f << std::endl;
     }
 
     convergence_iterations.push_back(f);
     ++nb_iterations;
   }
 
   this->computePressures();
 
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 Real BemGigipol::computeEquilibrium2init(Real epsilon, Real mean_pressure,
                                          Surface<Real>& init) {
 
   this->computeSpectralInfluenceOverDisplacement();
   this->surface_t = 0.;
   this->surface_r = 0.;
   this->search_direction = 0.;
   this->pold = 0.;
   this->surface_displacements = 0.;
   Real delta = 0.;
   Real Gold = 1.;
   Real f = 1e300;
   Real fPrevious = 1e300;
 
   UInt n = surface.size();
   UInt size = n * n;
 
   this->true_displacements = init;
 
   this->computeGaps();
 
 #pragma omp parallel for
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) < 0) {
       this->true_displacements(i) = this->surface(i);
     }
   }
   this->computeGaps();
 
   std::ofstream file("output.txt");
 
   convergence_iterations.clear();
   nb_iterations = 0;
   max_iterations = 200000;
 
   while (f > epsilon && nb_iterations < max_iterations) {
     fPrevious = f;
 
     this->functional->computeGradFU();
     this->search_direction = this->functional->getGradF();
 
     Real gbar = this->computeMeanPressuresInNonContact();
 
     this->search_direction += (2 * mean_pressure + gbar);
 
     Real G = this->computeG();
 
     this->updateT(G, Gold, delta);
     Real tau = this->computeTau(mean_pressure);
     tau = 0.01 * tau;
     this->old_displacements = this->true_displacements;
     Gold = G;
 #pragma omp parallel for
     for (unsigned int i = 0; i < size; ++i) {
 
       this->true_displacements(i) -= tau * this->surface_t(i);
     }
 
     // Projection on admissible space
     this->computeGaps();
 
 #pragma omp parallel for
     for (unsigned int i = 0; i < size; ++i) {
       if (this->gap(i) < 0) {
         this->true_displacements(i) = this->surface(i);
       }
     }
     this->computeGaps();
     delta = this->updateDisplacements(tau, 0.);
 
     // Projection on admissible space
     this->computeGaps();
 
     this->computePressures();
     f = this->computeStoppingCriterion();
 
     this->computeTractionsFromDisplacements();
     this->functional->computeGradFU();
     const Surface<Real>& gradF = this->functional->getGradF();
     Real min = SurfaceStatistics::computeMinimum(gradF);
 
     if (f < epsilon && epsilon < std::abs(min + mean_pressure))
       f = 3.;
 
     if (nb_iterations % dump_freq == 0) {
 
       std::cout << "G vaut " << G << std::endl;
       std::cout << "f vaut " << f << std::endl;
       std::cout << std::scientific << std::setprecision(10) << nb_iterations
                 << " " << f << " " << f - fPrevious << " " << 'G' << " " << G
                 << " " << std::endl;
 
       UInt n = surface.size();
       UInt size = n * n;
       Real crit = 0;
       Real disp_norm = 0;
       for (UInt i = 0; i < size; ++i) {
 
         crit += this->search_direction(i) * this->search_direction(i);
         disp_norm += (true_displacements(i) * true_displacements(i));
       }
 
       file << std::scientific << std::setprecision(10) << nb_iterations << " "
            << crit / disp_norm << " " << f << std::endl;
     }
 
     convergence_iterations.push_back(f);
     ++nb_iterations;
   }
 
   this->computePressures();
 
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::computeEquilibriuminit(Real epsilon, Real mean_displacement,
                                         Surface<Real>& init) {
 
   this->computeSpectralInfluenceOverDisplacement();
   this->surface_t = 0.;
   this->surface_r = 0.;
   this->search_direction = 0.;
   this->pold = 0.;
   this->surface_displacements = 0.;
   Real delta = 0.;
   Real Gold = 1.;
   Real f = 1e300;
   Real fPrevious = 1e300;
   UInt n = surface.size();
   UInt size = n * n;
 
   this->true_displacements = init;
   Real current_disp =
       SurfaceStatistics::computeAverage(this->true_displacements, 0);
   std::cout << "current disp " << current_disp << std::endl;
   if (current_disp == 0.) {
     true_displacements = mean_displacement;
     std::cout << "je reinitialise " << current_disp << std::endl;
   }
 
   this->computeGaps();
   this->optimizeToMeanDisplacement(mean_displacement);
 
   std::ofstream file("output.txt");
   this->computeGaps();
 
   convergence_iterations.clear();
   nb_iterations = 0;
   max_iterations = 100000;
 
   while (f > epsilon && nb_iterations < max_iterations) {
     fPrevious = f;
 
     this->functional->computeGradFU();
     this->search_direction = this->functional->getGradF();
     Real gbar = this->computeMeanPressuresInNonContact();
 
     this->search_direction -= gbar;
     Real G = this->computeG();
     this->updateT(G, Gold, delta);
     Real tau = this->computeTau();
 
     this->old_displacements = this->true_displacements;
     Gold = G;
 #pragma omp parallel for
     for (unsigned int i = 0; i < size; ++i) {
 
       this->true_displacements(i) -= tau * this->surface_t(i);
     }
 
     // Projection on admissible space
     this->computeGaps();
 
 #pragma omp parallel for
     for (unsigned int i = 0; i < size; ++i) {
       if (this->gap(i) < 0) {
         this->true_displacements(i) = this->surface(i);
       }
     }
     this->computeGaps();
     delta = this->updateDisplacements(tau, mean_displacement);
 
     // Projection on admissible space
     this->computeGaps();
     this->enforceMeanDisplacement(mean_displacement);
 
     this->computeGaps();
     this->computePressures();
     f = this->computeStoppingCriterion();
 
     if (nb_iterations % dump_freq == 0) {
 
       std::cout << "G vaut " << G << std::endl;
       std::cout << "f vaut " << f << std::endl;
       std::cout << std::scientific << std::setprecision(10) << nb_iterations
                 << " " << f << " " << f - fPrevious << " " << 'G' << " " << G
                 << " " << std::endl;
 
       UInt n = surface.size();
       UInt size = n * n;
       Real crit = 0;
       Real disp_norm = 0;
       for (UInt i = 0; i < size; ++i) {
 
         crit += this->search_direction(i) * this->search_direction(i);
         disp_norm += (true_displacements(i) * true_displacements(i));
       }
 
       file << std::scientific << std::setprecision(10) << nb_iterations << " "
            << crit / disp_norm << " " << f << std::endl;
     }
 
     convergence_iterations.push_back(f);
     ++nb_iterations;
   }
 
   this->computePressures();
 
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::computeStoppingCriterion() {
 
   UInt n = surface.size();
   UInt size = n * n;
 
   Real res = 0;
   Real t_sum =
       std::abs(SurfaceStatistics::computeSum(this->true_displacements));
   this->computeTractionsFromDisplacements();
   this->functional->computeGradFU();
   const Surface<Real>& gradF = this->functional->getGradF();
   Real min = SurfaceStatistics::computeMinimum(gradF);
 
 #pragma omp parallel for reduction(+ : res)
   for (UInt i = 0; i < size; ++i) {
     res +=
         std::abs((gradF(i) - min) * (this->true_displacements(i) - surface(i)));
 
     // res +=
     //   this->surface_tractions[i].real()
     //   *(surface_displacements[i].real() - surface[i].real());
   }
 
   return res / std::abs(t_sum);
 }
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::computeG() {
   unsigned int n = surface.size();
   unsigned int size = n * n;
   Real res = 0.;
 #pragma omp parallel for reduction(+ : res)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) > 0.) {
 
       Real val = this->search_direction(i);
       res += val * val;
     }
   }
 
   return res;
 }
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 Real BemGigipol::computeMeanPressuresInNonContact() {
   unsigned int n = surface.size();
   unsigned int size = n * n;
   Real res = 0.;
   UInt nb_contact = 0;
 #pragma omp parallel for reduction(+ : nb_contact, res)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) == 0.)
       continue;
     ++nb_contact;
     res += this->search_direction(i);
   }
 
   res /= nb_contact;
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigipol::optimizeToMeanDisplacement(Real imposed_mean) {
   Real target_value =
       imposed_mean - SurfaceStatistics::computeAverage(surface, 0);
   UInt n = surface.size();
   UInt size = n * n;
 
   // Initial guesses for upper and lower bound
   Real step_min = -10;
   Real step_max = 10;
 
   // Gaps for upper and lower bound
   Real gap_min = positiveGapAverage(step_min);
   Real gap_max = positiveGapAverage(step_max);
 
   UInt max_expansion = 8;
 
   // Expand bounds if necessary
   for (UInt i = 0; gap_max < target_value && i < max_expansion;
        i++, step_max *= 10)
     gap_max = positiveGapAverage(step_max);
 
   for (UInt i = 0; gap_min > target_value && i < max_expansion;
        i++, step_min *= 10)
     gap_min = positiveGapAverage(step_min);
 
   Real g = 0.;
   Real epsilon = 1e-12;
   Real step = 0;
 
   while (fabs(g - target_value) > epsilon) {
     step = (step_min + step_max) / 2.;
     g = positiveGapAverage(step);
     if (g > target_value)
       step_max = step;
     else if (g < target_value)
       step_min = step;
     else {
       step_max = step;
       step_min = step;
     }
   }
 
   step = (step_min + step_max) / 2.;
 
 #pragma omp parallel for
   for (UInt i = 0; i < size; i++) {
     gap(i) += step;
     if (gap(i) < 0)
       gap(i) = 0;
     true_displacements(i) = gap(i) + surface(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::positiveGapAverage(Real shift) {
   UInt n = surface.size();
 
   Real res = 0;
 #pragma omp parallel for reduction(+ : res)
   for (UInt i = 0; i < n * n; i++) {
     Real shifted_gap = gap(i) + shift;
     res += shifted_gap * (shifted_gap > 0);
   }
 
   return res / (n * n);
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigipol::updateT(Real G, Real Gold, Real delta) {
   unsigned int n = surface.size();
   unsigned int size = n * n;
   Real factor = delta * G / Gold;
 #pragma omp parallel for
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) == 0.)
       this->surface_t(i) = 0.;
     else {
       this->surface_t(i) *= factor;
       this->surface_t(i) += this->search_direction(i);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::computeTau() {
   const UInt n = surface.size();
   const UInt size = n * n;
 
   this->applyInverseInfluenceFunctions(surface_t, surface_r);
 
   Real rbar = 0;
   UInt nb_contact = 0;
 #pragma omp parallel for reduction(+ : nb_contact, rbar)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) == 0)
       continue;
     ++nb_contact;
     rbar += surface_r(i);
   }
   rbar /= nb_contact;
   surface_r -= rbar;
 
   Real tau_sum1 = 0., tau_sum2 = 0.;
 #pragma omp parallel for reduction(+ : tau_sum1, tau_sum2)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) == 0)
       continue;
     tau_sum1 += this->search_direction(i) * surface_t(i);
     tau_sum2 += surface_r(i) * surface_t(i);
   }
   Real tau = tau_sum1 / tau_sum2;
 
   return tau;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::computeTau(Real mean_pressure) {
   const UInt n = surface.size();
   const UInt size = n * n;
 
   this->applyInverseInfluenceFunctions(surface_t, surface_r);
 
   Real rbar = 0;
   UInt nb_contact = 0;
 #pragma omp parallel for reduction(+ : nb_contact, rbar)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) == 0)
       continue;
     ++nb_contact;
     rbar += surface_r(i);
   }
   rbar /= nb_contact;
   surface_r += 2 * mean_pressure + rbar;
 
   Real tau_sum1 = 0., tau_sum2 = 0.;
 #pragma omp parallel for reduction(+ : tau_sum1, tau_sum2)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) == 0)
       continue;
     tau_sum1 += this->search_direction(i) * surface_t(i);
     tau_sum2 += surface_r(i) * surface_t(i);
   }
   Real tau = tau_sum1 / tau_sum2;
 
   return tau;
 }
 
 /* -------------------------------------------------------------------------- */
 
-Real BemGigipol::updateDisplacements(Real tau, Real mean_displacement) {
+Real BemGigipol::updateDisplacements(Real tau, Real /*mean_displacement*/) {
   unsigned int n = surface.size();
   unsigned int size = n * n;
 
   // compute number of interpenetration without contact
   UInt nb_iol = 0;
 #pragma omp parallel for reduction(+ : nb_iol)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->search_direction(i) < 0 && this->gap(i) == 0.) {
 
       this->true_displacements(i) -= tau * this->search_direction(i);
       ++nb_iol;
     }
   }
 
   Real delta = 0;
   if (nb_iol > 0)
     delta = 0.;
   else
     delta = 1.;
 
   return delta;
 }
 
 /* -------------------------------------------------------------------------- */
 void BemGigipol::enforceMeanDisplacement(Real mean_displacement) {
   UInt n = surface.size();
   UInt size = n * n;
   Real moyenne_surface = SurfaceStatistics::computeAverage(this->surface, 0);
   Real average = SurfaceStatistics::computeAverage(this->true_displacements, 0);
   Real factor =
       (mean_displacement - moyenne_surface) / (average - moyenne_surface);
 
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     this->true_displacements(i) =
         factor * (this->true_displacements(i) - this->surface(i)) +
         this->surface(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigipol::computePressures() {
   this->computeTractionsFromDisplacements();
   this->functional->computeGradFU();
   const Surface<Real>& gradF = this->functional->getGradF();
   Real min = SurfaceStatistics::computeMinimum(gradF);
   this->surface_tractions -= min;
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
diff --git a/src/bem/bem_penalty.cpp b/src/bem/bem_penalty.cpp
index a21f133..820b425 100644
--- a/src/bem/bem_penalty.cpp
+++ b/src/bem/bem_penalty.cpp
@@ -1,200 +1,200 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 "bem_penalty.hh"
 #include "surface.hh"
 #include <fstream>
 #include <iomanip>
 #include <iostream>
 #include <sstream>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 void BemPenalty::computeTractionsFromDisplacements() {
   this->applyInverseInfluenceFunctions(this->true_displacements,
                                        this->surface_tractions);
 }
 
 Real BemPenalty::computeEquilibrium(Real epsilon, Real mean_displacement,
                                     Real penalization_parameter) {
 
   this->computeSpectralInfluenceOverDisplacement();
   this->search_direction = 0.;
   this->true_displacements = 0;
   this->true_displacements = mean_displacement;
   this->computeGaps();
   std::ofstream file("output.txt");
 
   Real f = 1e300;
 
   std::cout << "moyenne deplacement "
             << SurfaceStatistics::computeAverage(this->true_displacements, 0)
             << std::endl;
 
   convergence_iterations.clear();
   nb_iterations = 0;
 
   max_iterations = 5000;
 
   while (f > epsilon && nb_iterations++ < max_iterations) {
 
     this->functional->computeGradFU();
 
     this->computeSearchDirection(mean_displacement, penalization_parameter);
 
     Real alpha = 0.001;
 
     this->old_displacements = this->true_displacements;
 
     this->updateUnknown(alpha, mean_displacement);
 
     this->computeGaps();
 
     f = computeStoppingCriterion();
 
     if (nb_iterations % dump_freq == 0) {
       std::cout << std::scientific << std::setprecision(10) << nb_iterations
                 << " " << f << std::fixed << std::endl;
 
       this->computePressures(mean_displacement);
       Real orth = 0.;
       UInt n = surface.size();
       UInt size = n * n;
 
 #pragma omp parallel for reduction(+ : orth)
       for (UInt i = 0; i < size; ++i) {
 
         orth += std::abs(surface_tractions(i) *
                          (this->true_displacements(i) - surface(i)));
       }
 
       file << std::scientific << std::setprecision(10) << nb_iterations << " "
            << f << " " << orth << std::endl;
     }
 
     convergence_iterations.push_back(f);
   }
 
   this->computePressures(mean_displacement);
 
   file.close();
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemPenalty::computeStoppingCriterion() {
 
   Real crit = 0.;
   Real disp_norm = 0.;
 
   UInt n = surface.size();
   UInt size = n * n;
 
 #pragma omp parallel for reduction(+ : crit, disp_norm)
   for (UInt i = 0; i < size; ++i) {
 
     crit += this->search_direction(i) * this->search_direction(i);
     disp_norm += (true_displacements(i) * true_displacements(i));
   }
 
   return crit / disp_norm;
 }
 
 /* -------------------------------------------------------------------------- */
 
-void BemPenalty::computeSearchDirection(Real mean_displacement,
+void BemPenalty::computeSearchDirection(Real /*mean_displacement*/,
                                         Real penalization_parameter) {
   UInt n = surface.size();
   UInt size = n * n;
 
   const Surface<Real>& gradF = this->functional->getGradF();
 
 #pragma omp parallel for
   for (UInt i = 1; i < size; ++i) {
     this->search_direction(i) = gradF(i);
 
     if (gap(i) < 0) {
       this->search_direction(i) += penalization_parameter * gap(i);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemPenalty::computeOptimalStep() {
   this->applyInverseInfluenceFunctions(search_direction, surface_r);
   UInt n = surface.size();
   UInt size = n * n;
   Real numerator = 0., denominator = 0.;
 
 #pragma omp parallel for reduction(+ : numerator, denominator)
   for (UInt i = 0; i < size; ++i) {
     numerator += search_direction(i) * search_direction(i);
     denominator += surface_r(i) * search_direction(i);
   }
 
   Real alpha = numerator / denominator;
 
   return alpha;
 }
 
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 
 void BemPenalty::updateUnknown(Real alpha, Real mean_displacement) {
   UInt n = surface.size();
   UInt size = n * n;
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     this->true_displacements(i) -= alpha * this->search_direction(i);
   }
 
   Real moyenne = SurfaceStatistics::computeAverage(this->true_displacements, 0);
 
   for (UInt i = 0; i < size; ++i) {
     this->true_displacements(i) =
         this->true_displacements(i) - moyenne + mean_displacement;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
-void BemPenalty::computePressures(Real mean_displacement) {
+void BemPenalty::computePressures(Real /*mean_displacement*/) {
 
   this->computeTractionsFromDisplacements();
   this->functional->computeGradFU();
   // const Surface<Real> & gradF = this->functional->getGradF();
   // Real min = SurfaceStatistics::computeMinimum(gradF);
   this->surface_tractions -= this->surface_tractions(0);
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
diff --git a/src/bem/bem_uzawa.cpp b/src/bem/bem_uzawa.cpp
index cd4e67d..38c7b27 100644
--- a/src/bem/bem_uzawa.cpp
+++ b/src/bem/bem_uzawa.cpp
@@ -1,255 +1,255 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 "bem_uzawa.hh"
 #include "surface.hh"
 #include <fstream>
 #include <iomanip>
 #include <iostream>
 #include <sstream>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 void BemUzawa::computeTractionsFromDisplacements() {
   this->applyInverseInfluenceFunctions(this->true_displacements,
                                        this->surface_tractions);
 }
 
 Real BemUzawa::computeEquilibrium(Real epsilon, Real mean_displacement,
                                   Real penalization_parameter) {
 
   this->computeSpectralInfluenceOverDisplacement();
 
   this->search_direction = 0.;
   this->true_displacements = 0;
   this->true_displacements = mean_displacement;
   this->computeGaps();
   this->multipliers = 0;
 
   Real f_2 = 1e300;
   Real new_pen = 1e300;
   Real old_pen = 1e300;
 
   convergence_iterations.clear();
   UInt nb_iterations_1 = 0;
   UInt nb_iterations_2 = 0;
   std::ofstream file("output.txt");
 
   Real max_iterations_1 = 500;
   Real max_iterations_2 = 1000;
 
   while (new_pen > epsilon && nb_iterations_1++ < max_iterations_1) {
     std::cout << "iter ext " << nb_iterations_1 << std::endl;
 
     while (f_2 > 1e-12 && nb_iterations_2++ < max_iterations_2) {
       std::cout << "iter int " << nb_iterations_2 << std::endl;
       this->functional->computeGradFU();
       this->computeSearchDirection(mean_displacement, penalization_parameter);
 
       Real alpha = 1 / (10 * penalization_parameter);
       this->old_displacements = this->true_displacements;
       this->updateUnknown(alpha, mean_displacement);
       this->computeGaps();
 
       f_2 = computeStoppingCriterion();
 
       if (nb_iterations_2 % dump_freq == 0) {
         std::cout << std::scientific << std::setprecision(10) << nb_iterations
                   << " " << f_2 << std::fixed << std::endl;
 
         this->computePressures();
         Real orth = 0.;
         UInt n = surface.size();
         UInt size = n * n;
 
 #pragma omp parallel for reduction(+ : orth)
         for (UInt i = 0; i < size; ++i) {
 
           orth += std::abs(surface_tractions(i) *
                            (this->true_displacements(i) - surface(i)));
         }
 
         file << std::scientific << std::setprecision(10) << nb_iterations_2
              << " " << f_2 << " " << orth << std::endl;
       }
     }
 
     this->updateMultipliers(penalization_parameter);
 
     old_pen = computeInterpenetration(this->old_displacements);
     new_pen = computeInterpenetration(this->true_displacements);
     std::cout << "old penetration is  " << old_pen << std::endl;
     std::cout << "new penetration is  " << new_pen << std::endl;
 
     penalization_parameter =
         this->updatePenalization(penalization_parameter, old_pen, new_pen);
     // to avoid ill conditioned system
     if (penalization_parameter > 1.0e8)
       new_pen = 0.;
     nb_iterations_2 = 0;
     f_2 = 1e300;
     std::cout << "penalization is  " << penalization_parameter << std::endl;
   }
   this->computeGaps();
   this->computePressures();
 
   return new_pen;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemUzawa::computeStoppingCriterion() {
 
   Real crit = 0.;
   Real disp_norm = 0.;
 
   UInt n = surface.size();
   UInt size = n * n;
 
 #pragma omp parallel for reduction(+ : crit, disp_norm)
   for (UInt i = 0; i < size; ++i) {
 
     crit += (this->search_direction(i)) * (this->search_direction(i));
     disp_norm += (true_displacements(i) * true_displacements(i));
   }
 
   return crit / disp_norm;
 }
 
 /* -------------------------------------------------------------------------- */
 
-void BemUzawa::computeSearchDirection(Real mean_displacement,
+void BemUzawa::computeSearchDirection(Real /*mean_displacement*/,
                                       Real penalization_parameter) {
   UInt n = surface.size();
   UInt size = n * n;
 
   const Surface<Real>& gradF = this->functional->getGradF();
 
 #pragma omp parallel for
   for (UInt i = 1; i < size; ++i) {
     this->search_direction(i) = gradF(i);
 
     if (this->multipliers(i) + penalization_parameter * gap(i) <= 0) {
 
       this->search_direction(i) = this->search_direction(i) +
                                   this->multipliers(i) +
                                   penalization_parameter * gap(i);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemUzawa::computeOptimalStep() {
   this->applyInverseInfluenceFunctions(search_direction, surface_r);
   UInt n = surface.size();
   UInt size = n * n;
   Real numerator = 0., denominator = 0.;
 
 #pragma omp parallel for reduction(+ : numerator, denominator)
   for (UInt i = 0; i < size; ++i) {
     numerator += search_direction(i) * search_direction(i);
     denominator += surface_r(i) * search_direction(i);
   }
 
   Real alpha = numerator / denominator;
 
   return alpha;
 }
 
 /* -------------------------------------------------------------------------- */
 void BemUzawa::updateMultipliers(Real penalization_parameter) {
   UInt n = surface.size();
   UInt size = n * n;
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     if (this->multipliers(i) + penalization_parameter * gap(i) > 0.) {
       this->multipliers(i) = 0;
     } else {
       this->multipliers(i) += penalization_parameter * gap(i);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 Real BemUzawa::updatePenalization(Real penalization_parameter, Real old_pen,
                                   Real new_pen) {
   if (new_pen > 0.25 * old_pen) {
     penalization_parameter *= 5;
   }
 
   return penalization_parameter;
 }
 
 /* -------------------------------------------------------------------------- */
-Real BemUzawa::computeInterpenetration(Surface<Real>& displacements) {
+Real BemUzawa::computeInterpenetration(Surface<Real>& /*displacements*/) {
   UInt n = surface.size();
   UInt size = n * n;
   Real res = 0.;
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     if (0 > this->true_displacements(i) - surface(i)) {
       res += (this->true_displacements(i) - surface(i)) *
              (this->true_displacements(i) - surface(i));
     }
   }
 
   return sqrt(res);
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemUzawa::updateUnknown(Real alpha, Real mean_displacement) {
   UInt n = surface.size();
   UInt size = n * n;
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     this->true_displacements(i) -= alpha * this->search_direction(i);
   }
 
   Real moyenne = SurfaceStatistics::computeAverage(this->true_displacements, 0);
 
   for (UInt i = 0; i < size; ++i) {
     this->true_displacements(i) =
         this->true_displacements(i) - moyenne + mean_displacement;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemUzawa::computePressures() {
   this->computeTractionsFromDisplacements();
   this->functional->computeGradFU();
   // const Surface<Real> & gradF = this->functional->getGradF();
   // Real min = SurfaceStatistics::computeMinimum(gradF);
   this->surface_tractions -= this->surface_tractions(0);
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
diff --git a/src/core/grid.cpp b/src/core/grid.cpp
index e18317e..898f8d4 100644
--- a/src/core/grid.cpp
+++ b/src/core/grid.cpp
@@ -1,125 +1,126 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 "grid.hh"
 #include "tamaas.hh"
 #include <algorithm>
 #include <complex>
 #include <cstring>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid() : GridBase<T>() {
   this->n.fill(0);
   this->strides.fill(1);
   this->nb_components = 1;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid(const std::array<UInt, dim>& n, UInt nb_components,
                    value_type* data)
     : GridBase<T>() {
   this->n = n;
+  this->nb_components = nb_components;
   UInt size = this->computeSize();
   this->data.wrapMemory(data, size);
   this->computeStrides();
 }
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid(const Grid<T, dim>& o)
     : GridBase<T>(o), n(o.n), strides(o.strides) {}
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid(Grid<T, dim>&& o)
     : GridBase<T>(o), n(std::move(o.n)), strides(std::move(o.strides)) {}
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::resize(const std::array<UInt, dim>& n) {
   this->resize(n.begin(), n.end());
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::resize(const std::vector<UInt>& n) {
   TAMAAS_ASSERT(n.size() == dim, "Shape vector not matching grid dimensions");
   this->resize(n.begin(), n.end());
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::resize(std::initializer_list<UInt> n) {
   TAMAAS_ASSERT(n.size() == dim,
                 "Shape initializer list not matching grid dimensions");
   this->resize(std::begin(n), std::end(n));
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::computeStrides() {
   std::copy(n.begin() + 1, n.end(), strides.rbegin() + 2);
   strides[dim] = 1;
   strides[dim - 1] = this->nb_components;
   std::partial_sum(strides.rbegin(), strides.rend(), strides.rbegin(),
                    std::multiplies<UInt>());
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::printself(std::ostream& str) const {
   str << "Grid(" << dim << ", " << this->nb_components << ") {";
   for (auto& val : *this) {
     str << val << ", ";
   }
   str << "\b\b}";
 }
 
 /* -------------------------------------------------------------------------- */
 
 /// Class instanciation
 #define GRID_INSTANCIATE_TYPE(type)                                            \
   template class Grid<type, 1>;                                                \
   template class Grid<type, 2>;                                                \
   template class Grid<type, 3>
 
 GRID_INSTANCIATE_TYPE(Real);
 GRID_INSTANCIATE_TYPE(UInt);
 GRID_INSTANCIATE_TYPE(Complex);
 GRID_INSTANCIATE_TYPE(int);
 GRID_INSTANCIATE_TYPE(bool);
 GRID_INSTANCIATE_TYPE(unsigned long);
 #undef GRID_INSTANCIATE_TYPE
 
 __END_TAMAAS__
diff --git a/src/core/grid_view.hh b/src/core/grid_view.hh
index bbbaaba..cd888fb 100644
--- a/src/core/grid_view.hh
+++ b/src/core/grid_view.hh
@@ -1,161 +1,161 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 __GRID_VIEW_HH__
 #define __GRID_VIEW_HH__
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "tamaas.hh"
 #include <vector>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 /**
  * @brief View type on grid
  * This is a view on a *contiguous* chunk of data defined by a grid
  */
 template <template <typename, UInt> class Base, typename T, UInt base_dim,
 	  UInt dim>
 class GridView : public Base<T, dim> {
 public:
   using iterator = typename Base<T, dim>::iterator;
   using const_iterator = typename Base<T, dim>::const_iterator;
 
 public:
   /// Constructor
   GridView(GridBase<typename Base<T, dim>::value_type>& grid_base,
            const std::vector<UInt>& multi_index, Int component = -1);
   /// Move constructor
   GridView(GridView&& o)
       : Base<T, dim>(std::forward<Base<T, dim>>(o)),
 	grid(std::exchange(o.grid, nullptr)) {}
 
   /// Destructor
   ~GridView() override = default;
 
   UInt dataSize() const override { return this->computeSize(); }
 
   /// Iterators
-  iterator begin(UInt n = 1) override {
+  iterator begin(UInt /*n*/ = 1) override {
     return iterator(this->getInternalData(), this->strides.back());
   }
 
-  iterator end(UInt n = 1) override {
+  iterator end(UInt /*n*/ = 1) override {
     return iterator(this->getInternalData() +
                         this->dataSize() * this->strides.back(),
                     this->strides.back());
   }
 
-  const_iterator begin(UInt n = 1) const override {
+  const_iterator begin(UInt /*n*/ = 1) const override {
     return const_iterator(this->getInternalData(), this->strides.back());
   }
 
-  const_iterator end(UInt n = 1) const override {
+  const_iterator end(UInt /*n*/ = 1) const override {
     return const_iterator(this->getInternalData() +
                               this->dataSize() * this->strides.back(),
                           this->strides.back());
   }
 
 protected:
   Base<T, base_dim>* grid;
 };
 
 /* -------------------------------------------------------------------------- */
 template <template <typename, UInt> class Base, typename T, UInt base_dim,
 	  typename... Args>
 GridView<Base, T, base_dim, base_dim - sizeof...(Args)>
 make_view(Base<T, base_dim>& base, Args... indices) {
   std::vector<UInt> multi_index = {static_cast<UInt>(indices)...};
   return GridView<Base, T, base_dim, base_dim - sizeof...(Args)>(base,
 								 multi_index);
 }
 
 template <template <typename, UInt> class Base, typename T, UInt base_dim>
 GridView<Base, T, base_dim, base_dim>
 make_component_view(Base<T, base_dim>& base, UInt component) {
   std::vector<UInt> multi_index{};
   return GridView<Base, T, base_dim, base_dim>(base, multi_index, component);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template implementation */
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base, typename T, UInt base_dim,
           UInt dim>
 GridView<Base, T, base_dim, dim>::GridView(
     GridBase<typename Base<T, dim>::value_type>& grid_base,
     const std::vector<UInt>& multi_index, Int component)
     : Base<T, dim>(), grid(nullptr) {
   // static_assert(base_dim >= dim,
   // 		"view dimension must be >= than the base class");
 
   // Checking view type
   grid = dynamic_cast<Base<T, base_dim>*>(&grid_base);
 
   if (!grid)
     TAMAAS_EXCEPTION("given base type is incompatible with view");
 
   constexpr int dim_offset = base_dim - dim;
 
   if (dim_offset >= 0 and dim_offset != static_cast<Int>(multi_index.size()))
     TAMAAS_EXCEPTION("Number of blocked indices ("
                      << multi_index.size()
                      << ") does not match view dimension (" << base_dim
                      << " -> " << dim << ")");
 
   std::copy(grid->sizes().begin() + std::max(dim_offset, 0),
             grid->sizes().end(), this->n.begin() + std::max(-dim_offset, 0));
   std::copy(grid->getStrides().begin() + std::max(dim_offset, 0),
             grid->getStrides().end(),
             this->strides.begin() + std::max(-dim_offset, 0));
 
   if (dim_offset < 0)
     std::fill_n(this->n.begin(), -dim_offset, 1);
 
   if (component == -1) {  // view all components
     this->nb_components = grid->getNbComponents();
     component = 0;
   } else if (component >= 0 and
 	     component < static_cast<Int>(
 			     grid->getNbComponents())) {  // view one component
     this->nb_components = 1;
     this->strides.back() = grid->getNbComponents();
   } else
     TAMAAS_EXCEPTION("Components out of bounds for view");
 
   auto offset = std::inner_product(multi_index.begin(), multi_index.end(),
 				   grid->getStrides().begin(), component);
 
   this->data.wrapMemory(grid->getInternalData() + offset, this->computeSize());
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
 
 #endif  // __GRID_VIEW_HH__
diff --git a/src/core/loop.hh b/src/core/loop.hh
index 4e37c1b..9adfb46 100644
--- a/src/core/loop.hh
+++ b/src/core/loop.hh
@@ -1,226 +1,226 @@
 /**
  * @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 __LOOP_HH__
 #define __LOOP_HH__
 /* -------------------------------------------------------------------------- */
 #include "loops/apply.hh"
 #include "loops/loop_utils.hh"
 #include "tamaas.hh"
 #include <thrust/execution_policy.h>
 #include <thrust/for_each.h>
 #include <thrust/iterator/counting_iterator.h>
 #include <thrust/iterator/zip_iterator.h>
 #include <thrust/transform_reduce.h>
 #include <thrust/tuple.h>
 #include <type_traits>
 
 __BEGIN_TAMAAS__
 
 template <typename T>
 struct is_policy : std::false_type {};
 template <>
 struct is_policy<thrust::detail::host_t> : std::true_type {};
 template <>
 struct is_policy<thrust::detail::device_t> : std::true_type {};
 template <>
 struct is_policy<const thrust::detail::host_t> : std::true_type {};
 template <>
 struct is_policy<const thrust::detail::device_t> : std::true_type {};
 template <>
 struct is_policy<const thrust::detail::host_t&> : std::true_type {};
 template <>
 struct is_policy<const thrust::detail::device_t&> : std::true_type {};
 
 /**
  * @brief Singleton class for automated loops using lambdas
 
  * This class is sweet candy :) It provides abstraction of the paralelism
  * paradigm used in loops and allows simple and less erro-prone loop syntax,
  * with minimum boiler plate. I love it <3
  */
 class Loop {
 public:
   /// Backends enumeration
   enum backend {
     omp,   ///< [OpenMP](http://www.openmp.org/specifications/) backend
     cuda,  ///< [Cuda](http://docs.nvidia.com/cuda/index.html) backend
   };
 
   /// Helper class to count iterations within lambda-loop
   template <typename T>
   class arange {
   public:
     using it_type = thrust::counting_iterator<T>;
     arange(T start, T size) : start(start), range_size(size) {}
-    it_type begin(UInt i = 0) const { return it_type(T(start)); }
-    it_type end(UInt i = 0) const { return it_type(range_size); }
+    it_type begin(UInt /*i*/ = 1) const { return it_type(T(start)); }
+    it_type end(UInt /*i*/ = 1) const { return it_type(range_size); }
     UInt getNbComponents() const { return 1; }
 
   private:
     T start, range_size;
   };
 
   template <typename T>
   static arange<T> range(T size) {
     return arange<T>(0, size);
   }
 
   template <typename T, typename U>
   static arange<T> range(U start, T size) {
     return arange<T>(start, size);
   }
 
   /// Loop functor over any number of grids
   template <typename Functor, typename... Grids>
   static void loop(Functor&& func, Grids&&... containers);
 
   /// Strided loop over any number of grids with parallel policy
   template <typename DerivedPolicy, typename Functor, typename... Grids>
   static void stridedLoop(const thrust::execution_policy<DerivedPolicy>& policy,
                           Functor&& func, Grids&&... containers) {
     stridedLoopImpl(policy, std::forward<Functor>(func),
                     std::forward<Grids>(containers)...);
   }
 
   /// Strided loop over any number of grids
   template <typename Functor, typename... Grids>
   static typename std::enable_if<not is_policy<Functor>::value, void>::type
   stridedLoop(Functor&& func, Grids&&... containers) {
     stridedLoopImpl(thrust::device, std::forward<Functor>(func),
                     std::forward<Grids>(containers)...);
   }
 
 private:
   /// Implementation of strided loop overloads
   template <typename DerivedPolicy, typename Functor, typename... Grids>
   static void
   stridedLoopImpl(const thrust::execution_policy<DerivedPolicy>& policy,
                   Functor&& func, Grids&&... containers);
 
 public:
   /// Reduce over any number of grids
   template <operation op, typename Functor, typename... Grids>
   static auto reduce(Functor&& func, Grids&&... containers)
       -> decltype(func(containers(0)...));
 
   /// Strided reduce over any number of grids
   template <operation op, typename Functor, typename... Grids>
   static auto stridedReduce(Functor&& func, Grids&&... containers)
       -> decltype(func(containers(0)...));
 
   /// Constructor
   Loop() = delete;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Template implementation                                                    */
 /* -------------------------------------------------------------------------- */
 
 template <typename Functor, typename... Grids>
 void Loop::loop(Functor&& func, Grids&&... containers) {
   auto begin =
       thrust::make_zip_iterator(thrust::make_tuple(containers.begin()...));
   auto end = thrust::make_zip_iterator(thrust::make_tuple(containers.end()...));
 
   thrust::for_each(begin, end, detail::ApplyFunctor<Functor>(func));
 
 #ifdef USE_CUDA
   cudaDeviceSynchronize();
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename DevicePolicy, typename Functor, typename... Grids>
 void Loop::stridedLoopImpl(
     const thrust::execution_policy<DevicePolicy>& /*policy*/, Functor&& func,
     Grids&&... containers) {
   auto begin = thrust::make_zip_iterator(
       thrust::make_tuple(containers.begin(containers.getNbComponents())...));
   auto end = thrust::make_zip_iterator(
       thrust::make_tuple(containers.end(containers.getNbComponents())...));
 
   thrust::for_each(begin, end, detail::ApplyFunctor<Functor>(func));
 #ifdef USE_CUDA
   cudaDeviceSynchronize();
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <operation op, typename Functor, typename... Grids>
 auto Loop::reduce(Functor&& func, Grids&&... containers)
     -> decltype(func(containers(0)...)) {
   auto begin =
       thrust::make_zip_iterator(thrust::make_tuple(containers.begin()...));
   auto end = thrust::make_zip_iterator(thrust::make_tuple(containers.end()...));
 
   using reduce_type = decltype(func(containers(0)...));
   using apply_type = detail::ApplyFunctor<Functor, reduce_type>;
 
   auto red_helper = detail::reduction_helper<op, apply_type>(apply_type(func));
   auto result = thrust::reduce(
       begin, end, red_helper.template init<reduce_type>(), red_helper);
 
 #ifdef USE_CUDA
   cudaDeviceSynchronize();
 #endif
   return result;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <operation op, typename Functor, typename... Grids>
 auto Loop::stridedReduce(Functor&& func, Grids&&... containers)
     -> decltype(func(containers(0)...)) {
   auto begin = thrust::make_zip_iterator(
       thrust::make_tuple(containers.begin(containers.getNbComponents())...));
   auto end = thrust::make_zip_iterator(
       thrust::make_tuple(containers.end(containers.getNbComponents())...));
 
   using reduce_type = decltype(func(containers(0)...));
   using apply_type = detail::ApplyFunctor<Functor, reduce_type>;
 
   auto red_helper = detail::reduction_helper<op, apply_type>(apply_type(func));
 
   auto result = thrust::reduce(
       begin, end, red_helper.template init<reduce_type>(), red_helper);
 
 #ifdef USE_CUDA
   cudaDeviceSynchronize();
 #endif
   return result;
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
 
 #undef EXEC_CASE_MACRO
 #undef REDUCE_CASE_MACRO
 
 #endif  // __LOOP_HH__
diff --git a/src/core/surface.cpp b/src/core/surface.cpp
index c0038d4..f9b8ccf 100644
--- a/src/core/surface.cpp
+++ b/src/core/surface.cpp
@@ -1,427 +1,427 @@
 /**
  * @file
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 <fstream>
 #include <sstream>
 #ifdef USE_CUDA
 #include <thrust/random.h>
 #else
 #include <random>
 #endif
 /* -------------------------------------------------------------------------- */
 #include "surface.hh"
 #include "surface_statistics.hh"
 /* -------------------------------------------------------------------------- */
 #ifdef USING_IOHELPER
 #include "dumper_paraview.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 std::string fname = "surface_height.pdf";
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::expandByLinearInterpolation() {
   UInt n = this->size();
   Surface* old_s = new Surface(n, 1.);
   int new_n = n * 2;
   *old_s = *this;
   this->setGridSize(new_n);
   for (UInt i = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       (*this)(i, j) = (*old_s)(i, j);
     }
   }
   delete old_s;
 }
 
 Real sinc(Real x) { return sin(M_PI * x) / (M_PI * x); }
 /* -------------------------------------------------------------------------- */
 template <typename T>
-void Surface<T>::expandByShannonInterpolation(int factor) {
+void Surface<T>::expandByShannonInterpolation(int /*factor*/) {
 
   /* /!\ this->FFTTransform(); /!\ */
   TAMAAS_EXCEPTION("Shannon interpolation needs old FFTTransform interface");
 #if 0
   UInt n = this->size();
   UInt N = std::pow(2,factor)*this->size();
   std::cout << N << " " << factor << std::endl;
   Surface s = Surface(N,1.);
 
 
   for (UInt i = 0; i < n/2; ++i)
     for (UInt j = 0; j < n/2; ++j) {
       s(i,j) = (*this)(i,j);
     }
 
   for (UInt i = n/2; i < n; ++i)
     for (UInt j = 0; j < n/2; ++j) {
       s(N-n+i,j) = (*this)(i,j);
     }
 
   for (UInt i = n/2; i < n; ++i)
     for (UInt j = n/2; j < n; ++j) {
       s(N-n+i,N-n+j) = (*this)(i,j);
     }
 
   for (UInt i = 0; i < n/2; ++i)
     for (UInt j = n/2; j < n; ++j) {
       s(i,N-n+j) = (*this)(i,j);
     }
 
   //  s.dumpToTextFile<false,true>("temp.csv");
 
   /* /!\ s.FFTITransform(); /!\ */
 
   this->setGridSize(N);
   Real f = std::pow(2,factor);
   f = f*f;
 
   for (UInt i = 0; i < N*N; ++i)
     (*this)(i) = f*s(i);
 
 
 //   int new_n = 2*n;
 //   Surface s(*this);
 //   this->setGridSize(new_n);
 
 //   // put the original signal with zeros for to be interpolated points
 //   for (int i=0; i<n; i++)
 //     for (int j=0; j<n; j++)
 //       (*this)(2*i,j*2) = s(i,j);
 
 //   // do shannon interpolation
 //   int rcut = 64;
 
 //   //interpolate lines
 //   for (int i=0; i<new_n; i++) {
 //     for (int j=0; j<new_n; j++) {
 //       if (i%2 == 0) continue;
 //       // Real loop on local neighbors to interpolate
 //       for (int k=-rcut; k<rcut+1; k++) {
 // 	int i2 = i + k;
 // 	if (i2 < 0) i2 += new_n;
 // 	if (i2 >= new_n) i2 -= new_n;
 // 	if (i2%2) continue;
 
 // 	Real d = fabs(k);
 // 	s[i+j*new_n] += s(i2,j)*sinc(d/2);
 //       }
 //     }
 //   }
 
 //   //interpolate columns
 //   for (int i=0; i<new_n; i++) {
 //     for (int j=0; j<new_n; j++) {
 //       if (j%2 == 0) continue;
 //       // Real loop on local neighbors to interpolate
 //       for (int k=-rcut; k<rcut+1; k++) {
 // 	int j2 = j + k;
 // 	if (j2 < 0) j2 += new_n;
 // 	if (j2 >= new_n) j2 -= new_n;
 // 	if (j2%2) continue;
 
 // 	Real d = fabs(k);
 // 	s[i+j*new_n] += complex(s[i+j2*new_n].real()*sinc(d/2),
 // 				s[i+j2*new_n].imag()*sinc(d/2));
 //       }
 //     }
 //   }
 
 
 // //   for (int i=0; i<n; i++) {
 // //     for (int j=0; j<n; j++) {
 // //       s[2*i+j*2*new_n].real() = 0;
 // //       s[2*i+j*2*new_n].im = 0;
 // //     }
 // //   }
 #endif
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByMirroring() {
   UInt n = this->size();
   UInt new_n = n * 2;
   Surface s(*this);
   this->setGridSize(new_n);
 
   for (UInt i = 0; i < new_n; i++) {
     for (UInt j = 0; j < new_n; j++) {
       UInt i1 = i;
       if (i1 > n)
         i1 = new_n - i1;
       UInt j1 = j;
       if (j1 > n)
         j1 = new_n - j1;
       if (i == n || j == n) {
         (*this)(i, j) = 0.;
         continue;
       }
       (*this)(i, j) = s(i1, j1);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByPuttingZeros() {
   UInt n = this->size();
   UInt new_n = n * 2;
   Surface s(*this);
   this->setGridSize(new_n);
 
   for (UInt i = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       (*this)(2 * i, j * 2) = s(i, j);
     }
   }
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::contractByTakingSubRegion() {
   UInt n = this->size();
   Surface s(*this);
   this->setGridSize(n / 2);
 
   for (UInt i = 0; i < n / 2; i++) {
     for (UInt j = 0; j < n / 2; j++) {
       (*this)(i, j) = s(i, j);
     }
   }
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::contractByIgnoringMidPoints() {
   UInt n = this->size();
   Surface s(*this);
   this->setGridSize(n / 2);
 
   for (UInt i = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       if (i % 2)
         continue;
       if (j % 2)
         continue;
       (*this)(i / 2, j / 2) = s(i, j);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 // void Surface<T>::dumpToTextFile(string filename, bool pr) {
 
 //   cout << "Writing surface file " << filename << endl;
 //   ofstream file(filename.c_str());
 //   file.precision(15);
 
 //   if (pr)
 //   {
 //     for (int i=0; i<n; i++) {
 //       for (int j=0; j<n; j++) {
 //         file << i << " " << j << " ";
 //         if (heights[i+j*n].real() < 1e-10)
 //           file << scientific << -1e10 << " " << 0;
 //         else
 //           file << scientific << heights[i+j*n].real() << " " <<
 //           heights[i+j*n].im;
 //         file << endl;
 //       }
 //       file << endl;
 //     }
 //   }
 //   else
 //   {
 //     for (int i=0; i<n; i++) {
 //       for (int j=0; j<n; j++) {
 //         file << i << " " << j << " ";
 //         file << scientific << heights[i+j*n].real() << " " <<
 //         heights[i+j*n].im; file << endl;
 //       }
 //     }
 //   }
 //   file.close();
 
 // }
 
 /* -------------------------------------------------------------------------- */
 
 // Surface * Surface<T>::loadFromTextFile(std::string filename) {
 
 // std::string line;
 // std::ifstream myfile (filename.c_str());
 
 // if (!myfile.is_open()){
 //   SURFACE_FATAL("cannot open file " << filename);
 // }
 
 // Surface * my_surface = nullptr;
 // UInt n = 0;
 // if (myfile.is_open())
 //   {
 //     while (! myfile.eof() )
 // 	{
 // 	  getline (myfile,line);
 // 	  if (line[0] == '#')
 // 	    continue;
 // 	  if (line.size() == 0)
 // 	    continue;
 
 // 	  std::stringstream s;
 // 	  s << line;
 // 	  int i,j;
 // 	  s >> i;
 // 	  s >> j;
 // 	  if (n < i) n = i;
 // 	  if (n < j) n = j;
 // 	}
 //     ++n;
 //     std::cout << "read " << n << " x " << n << " grid points" << std::endl;
 
 //     myfile.clear();              // forget we hit the end of file
 //     myfile.seekg(0, std::ios::beg);   // move to the start of the file
 
 //     my_surface = new Surface(n);
 
 //     while (! myfile.eof() )
 // 	{
 // 	  getline (myfile,line);
 // 	  if (line[0] == '#')
 // 	    continue;
 // 	  if (line.size() == 0)
 // 	    continue;
 
 // 	  std::stringstream s;
 // 	  s << line;
 // 	  int i,j;
 // 	  Real h,hIm;
 // 	  s >> i;
 // 	  s >> j;
 // 	  s >> h;
 // 	  s >> hIm;
 // 	  my_surface->heights[i+j*n] = complex(h,hIm);
 // 	}
 //     myfile.close();
 //   }
 //   return my_surface;
 // }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
-void Surface<T>::dumpToParaview(std::string filename) const {
+void Surface<T>::dumpToParaview(std::string /*filename*/) const {
 #ifdef USING_IOHELPER
 
   UInt n = this->size();
   const Surface& data = *this;
   /* build 3D coordinates array to be usable with iohelper */
   //  cerr << "allocate coords for paraview" << endl;
   Real* coords = new Real[n * n * 3];
   // Real *  coords = new Real[(2*n-1)*(2*n-1)*3];
   // cerr << "allocate zcoords for paraview" << endl;
   Real* zcoords = new Real[n * n * 3];
   // Real *  zcoords = new Real[(2*n-1)*(2*n-1)*3];
   // cerr << "allocate zcoords2 for paraview" << endl;
   Real* zcoords2 = new Real[n * n * 3];
   //  Real *  zcoords2 = new Real[(2*n-1)*(2*n-1)*3];
   int index = 0;
   // for (int i=-n+1; i<n; i++) {
   //   for (int j=-n+1; j<n; j++) {
   for (int i = 0; i < n; i++) {
     for (int j = 0; j < n; j++) {
       coords[3 * index] = 1. * i / n;
       coords[3 * index + 1] = 1. * j / n;
       coords[3 * index + 2] = 0;
       zcoords[3 * index] = 0;
       zcoords[3 * index + 1] = 0;
       zcoords2[3 * index] = 0;
       zcoords2[3 * index + 1] = 0;
       int i1 = (i % n);
       int j1 = (j % n);
       if (i1 < 0)
         i1 += n;
       if (j1 < 0)
         j1 += n;
       zcoords[3 * index + 2] = data(i1, j1).real();
       zcoords2[3 * index + 2] = data(i1, j1).imag();
       ++index;
     }
   }
 
   iohelper::DumperParaview dumper;
   dumper.setPoints(coords, 3, index, filename.c_str());
   dumper.setConnectivity(nullptr, iohelper::POINT_SET, 1, iohelper::FORTRAN_MODE);
   dumper.addNodeDataField("data.re", zcoords, 3, index);
   dumper.addNodeDataField("data.im", zcoords2, 3, index);
   dumper.setPrefix("./paraview");
   dumper.init();
   dumper.dump();
 
   delete[] coords;
   delete[] zcoords;
   delete[] zcoords2;
 #endif  // USING_IOHELPER
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::recenterHeights() {
   UInt n = this->size();
   Real zmean = SurfaceStatistics::computeAverage(*this);
 
   for (UInt i = 0; i < n * n; ++i)
     (*this)(i) = (*this)(i)-zmean;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::generateWhiteNoiseSurface(Surface<Real>& sources,
                                            long random_seed) {
   // set the random gaussian generator
   random_engine gen(random_seed);
   normal_distribution<Real> gaussian_generator;
 
   UInt n = sources.size();
 
   for (UInt i = 0; i < n; ++i)
     for (UInt j = 0; j < n; ++j) {
       sources(i, j) = gaussian_generator(gen);
     }
 }
 /* -------------------------------------------------------------------------- */
 
 template class Surface<Real>;
 template class Surface<unsigned long>;
 template class Surface<UInt>;
 template class Surface<int>;
 template class Surface<bool>;
 
 __END_TAMAAS__
diff --git a/src/core/surface_complex.hh b/src/core/surface_complex.hh
index 86013bf..176b69a 100644
--- a/src/core/surface_complex.hh
+++ b/src/core/surface_complex.hh
@@ -1,163 +1,163 @@
 /**
  * @file
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 SURFACE_COMPLEX_H
 #define SURFACE_COMPLEX_H
 /* -------------------------------------------------------------------------- */
 #include "grid_hermitian.hh"
 #include "surface.hh"
 #include <thrust/complex.h>
 
 __BEGIN_TAMAAS__
 
 /**
  * @deprecated should use GridHermitian instead
  * @brief Square 2D complex surface
  */
 template <typename T>
 class SurfaceComplex : public GridHermitian<T, 2> {
 
 public:
   SurfaceComplex(UInt a, Real L);
   using GridHermitian<T, 2>::GridHermitian;
   using GridHermitian<T, 2>::operator=;
   using GridHermitian<T, 2>::operator();
 
   //! make a surface real by taking norm as real part and imaginary 0
   void makeItRealBySquare();
   //! make it real by summing the imaginary and real parts
   void makeItRealBySum();
   //! make a surface real by taking norm as real part and imaginary 0
   void makeItRealByAbs();
   //! get real part
   Surface<T> real() const;
   //! get imaginary part
   Surface<T> imag() const;
 
   UInt size() const { return this->n[0]; }
   UInt getL() const { return 1.; }
 
   void setGridSize(UInt s);
 
   virtual const SurfaceComplex<T>& getWrappedNumpy() { return *this; };
 
 #define SURF_OPERATOR(op)                                                      \
   inline void operator op(const Surface<T>& other);                            \
   inline void operator op(const SurfaceComplex<T>& other)
 
   SURF_OPERATOR(+=);
   SURF_OPERATOR(*=);
   SURF_OPERATOR(-=);
   SURF_OPERATOR(/=);
 #undef SURF_OPERATOR
 };
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void SurfaceComplex<T>::setGridSize(UInt s) {
   this->n[0] = s;
   this->n[1] = s / 2 + 1;
   this->resize(this->n);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 Surface<T> SurfaceComplex<T>::real() const {
   Surface<T> res(this->size(), this->getL());
   UInt n = this->size();
   for (UInt i = 0; i < n; ++i)
     for (UInt j = 0; j < n; ++j) {
       res(i, j) = (*this)(i, j).real();
     }
   return res;
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 Surface<T> SurfaceComplex<T>::imag() const {
   Surface<T> res(this->size(), this->getL());
   UInt n = this->size();
   for (UInt i = 0; i < n; ++i)
     for (UInt j = 0; j < n; ++j) {
       res(i, j) = (*this)(i, j).imag();
     }
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void SurfaceComplex<T>::makeItRealBySquare() {
   Loop::loop([] CUDA_LAMBDA(complex<T> & x) { x *= thrust::conj(x); }, *this);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void SurfaceComplex<T>::makeItRealByAbs() {
   Loop::loop([] CUDA_LAMBDA(complex<T> & x) { x = thrust::abs(x); }, *this);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void SurfaceComplex<T>::makeItRealBySum() {
   Loop::loop([] CUDA_LAMBDA(complex<T> & x) { x = x.real() + x.imag(); },
              *this);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
-SurfaceComplex<T>::SurfaceComplex(UInt a, Real L) : GridHermitian<T, 2>() {
+SurfaceComplex<T>::SurfaceComplex(UInt a, Real /*L*/) : GridHermitian<T, 2>() {
   this->n = GridHermitian<T, 2>::hermitianDimensions({{a, a}});
   this->resize(this->n);
 }
 
 /* -------------------------------------------------------------------------- */
 #define SURF_OP_IMPL(op)                                                       \
   template <typename T>                                                        \
   inline void SurfaceComplex<T>::operator op(const Surface<T>& other) {        \
     TAMAAS_ASSERT(this->n[0] == other.sizes()[0] &&                            \
                       this->n[1] == other.sizes()[1] / 2 + 1,                  \
                   "surface size does not match");                              \
     _Pragma("omp parallel for") for (UInt i = 0; i < this->n[0];               \
                                      i++) for (UInt j = 0; j < this->n[1];     \
                                                j++) (*this)(i, j)op other(i,   \
                                                                           j);  \
   }                                                                            \
   template <typename T>                                                        \
   inline void SurfaceComplex<T>::operator op(const SurfaceComplex<T>& other) { \
     Grid<complex<T>, 2>::operator op(other);                                   \
   }
 
 SURF_OP_IMPL(+=);
 SURF_OP_IMPL(*=);
 SURF_OP_IMPL(-=);
 SURF_OP_IMPL(/=);
 #undef SURF_OP_IMPL
 
 __END_TAMAAS__
 
 #endif /* SURFACE_COMPLEX_H */
diff --git a/src/core/surface_statistics.hh b/src/core/surface_statistics.hh
index 0c1287b..92aa42d 100644
--- a/src/core/surface_statistics.hh
+++ b/src/core/surface_statistics.hh
@@ -1,721 +1,721 @@
 /**
  * @file
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 __SURFACE_STATISTICS_HH__
 #define __SURFACE_STATISTICS_HH__
 /* -------------------------------------------------------------------------- */
 #include "fft_plan_manager.hh"
 #include "fftransform.hh"
 #include "surface.hh"
 #include <thrust/complex.h>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /**
  * @deprecated should move code to Statistics class instead
  * @brief Statistics calculations on surfaces
  */
 class SurfaceStatistics {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   SurfaceStatistics(){};
   virtual ~SurfaceStatistics(){};
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   inline static Real computeMaximum(const Surface<Real>& s);
   inline static Real computeMinimum(const Surface<Real>& s);
-  inline static UInt computeAverage(const Surface<UInt>& s, UInt epsilon = 0) {
+  inline static UInt computeAverage(const Surface<UInt>& /*s*/, UInt /*epsilon*/) {
     SURFACE_FATAL("TO IMPLEMENT");
   };
 
-  inline static unsigned long computeAverage(const Surface<unsigned long>& s,
-                                             unsigned long epsilon = 0) {
+  inline static unsigned long computeAverage(const Surface<unsigned long>& /*s*/,
+                                             unsigned long /*epsilon*/) {
     SURFACE_FATAL("TO IMPLEMENT");
   };
 
   inline static Real
   computeAverage(const Surface<Real>& s,
                  Real epsilon = std::numeric_limits<Real>::lowest()) {
     return computeAverage<Real>(s, epsilon);
   }
 
   template <typename T>
   inline static T computeAverage(const Surface<T>& s,
                                  T epsilon = std::numeric_limits<T>::lowest());
   inline static Real
   computeStdev(const Grid<Real, 2>& s,
                Real epsilon = std::numeric_limits<Real>::lowest());
 
   inline static Real computeSkewness(const Surface<Real>& surface, Real avg,
                                      Real std);
   inline static Real computeKurtosis(const Surface<Real>& surface, Real avg,
                                      Real std);
   inline static Real computeRMSSlope(const Surface<Real>& surface);
   inline static Real computeSpectralRMSSlope(Surface<Real>& surface);
   inline static Real computeSpectralMeanCurvature(Surface<Real>& surface);
   inline static Real computeSpectralStdev(Surface<Real>& surface);
   //  inline static void computeStatistics(Surface<Real> & surface);
   inline static std::vector<Real> computeMoments(Surface<Real>& surface);
   inline static Real computeAnalyticFractalMoment(UInt order, UInt k1, UInt k2,
                                                   Real hurst, Real rms, Real L);
   //! compute perimeter
   inline static Real computePerimeter(const Surface<Real>& surface,
                                       Real eps = 1e-10);
   //! get real contact area
   inline static Real computeContactArea(const Surface<Real>& surface);
   //! get real contact area
   inline static Real computeContactAreaRatio(const Surface<Real>& surface);
 
   template <bool consider_zeros = false>
   inline static std::vector<int>
   computeDistribution(const Surface<Real>& surface, unsigned int n_bins,
                       Real min_value = std::numeric_limits<Real>::max(),
                       Real max_value = std::numeric_limits<Real>::lowest());
 
   inline static std::vector<Complex>
   computeSpectralDistribution(const Surface<Real>& surface, unsigned int n_bins,
                               UInt cutoff = std::numeric_limits<UInt>::max());
 
   //! compute average amplitude
   inline static Real computeAverageAmplitude(const Surface<Real>& surface);
   inline static Complex
   computeAverageAmplitude(const SurfaceComplex<Real>& surface);
 
   //! sum all heights
   inline static Real computeSum(const Surface<Real>& surface);
   //! compute auto correlation
   inline static Surface<Real> computeAutocorrelation(Surface<Real>& surface);
   //! compute powerspectrum
   inline static SurfaceComplex<Real>
   computePowerSpectrum(Surface<Real>& surface);
 };
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 
 //#include "surface_statistics_inline_impl.cc"
 
 /* -------------------------------------------------------------------------- */
 
 template <bool consider_zeros>
 std::vector<int>
 SurfaceStatistics::computeDistribution(const Surface<Real>& surface,
                                        unsigned int n_bins, Real min_value,
                                        Real max_value) {
 
   if (n_bins == 0)
     SURFACE_FATAL("badly set number of bins");
 
   Real min_heights;
   Real max_heights;
 
   std::vector<int> bins;
   bins.resize(n_bins);
   bins.assign(n_bins, 0);
 
   min_heights = min_value;
   max_heights = max_value;
 
   if (max_value == -1. * std::numeric_limits<Real>::max())
     max_heights = computeMaximum(surface);
   if (min_value == std::numeric_limits<Real>::max())
     min_heights = computeMinimum(surface);
 
   int counter = 0;
   const int n = surface.size();
   for (int i = 0; i < n * n; ++i) {
     Real value = surface(i);
     if (!consider_zeros && value < zero_threshold)
       continue;
     if (value > max_heights || value < min_heights)
       SURFACE_FATAL("out of bound value: "
                     << min_heights << " <? " << value << " <? " << max_heights
                     << " : bounds probably badly chosen");
     Real ind = (value - min_heights) / (max_heights - min_heights);
     unsigned int index = ind * n_bins;
     if (index == n_bins)
       --index;
     ++bins[index];
     ++counter;
   }
   if (consider_zeros && counter != n * n)
     SURFACE_FATAL("problem in counting points");
 
   //  Real offset = (max_heights - min_heights)/n_bins;
   // Real x = 0;
   //  Real f1 = 0, f2 = 0;
   //  for (unsigned int i = 0 ; i < n_bins ; ++i)
   //  {
   //    x  = i*offset+min_heights;
   // f1 = bins[i]/Real(n*n);
   // f2 = bins[i]/Real(counter);
   //}
   Real min_non_zero = std::numeric_limits<Real>::max();
 #pragma omp parallel for reduction(min : min_non_zero)
   for (int i = 0; i < n * n; ++i) {
     Real value = surface(i);
     if (value < min_non_zero && value > 1e-14)
       min_non_zero = value;
   }
   return bins;
 }
 /* -------------------------------------------------------------------------- */
 std::vector<Complex> SurfaceStatistics::computeSpectralDistribution(
     const Surface<Real>& surface, unsigned int n_bins, UInt cutoff) {
 
   if (n_bins == 0)
     SURFACE_FATAL("badly set number of bins");
 
   const int n = surface.size();
 
   Real min_value = SurfaceStatistics::computeMinimum(surface);
   Real max_value = SurfaceStatistics::computeMaximum(surface);
 
   std::vector<Complex> distrib_fft(n_bins);
   const int N = n_bins;
   Real L = max_value - min_value;
 
 #pragma omp parallel for
   for (int k = 0; k < N; ++k)
     distrib_fft[k] = 0.;
 
 #pragma omp parallel for
   for (int k = 1; k < N / 2; ++k) {
     for (int i = 0; i < n * n; ++i) {
       Real value = surface(i);
       Real angle = -2. * M_PI * k * (value - min_value) / L;
       Complex phase = std::polar(1., angle);
       distrib_fft[k] += phase;
       distrib_fft[N - k] += phase;
     }
   }
 
   for (int i = 0; i < n * n; ++i) {
     Real value = surface(i);
     Real angle = -2. * M_PI * N / 2 * (value - min_value) / L;
     Complex phase = std::polar(1., angle);
     distrib_fft[N / 2] += phase;
   }
 
 #pragma omp parallel for
   for (int k = 1; k < N; ++k)
     distrib_fft[k] /= n * n;
 
   distrib_fft[0] = std::polar(1., 0.);
 
 #pragma omp parallel for
   for (int k = 1; k < N / 2; ++k) {
     if (k >= (int)cutoff) {
       distrib_fft[k] = 0.;
       distrib_fft[N - k] = 0.;
     }
   }
 
     // proceed to inverse fourier transform
 
 #pragma omp parallel for
   for (int k = 0; k < N; ++k) {
     distrib_fft[k] = conj(distrib_fft[k]);
   }
 
   fftw_complex* in = (fftw_complex*)&distrib_fft[0];
   fftw_plan p = fftw_plan_dft_1d(N, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
   fftw_execute(p);
   fftw_destroy_plan(p);
 
   std::vector<Complex> distrib;
   distrib.resize(N);
 
 #pragma omp parallel for
   for (int k = 0; k < N; ++k) {
     distrib[k] = conj(distrib_fft[k]);
   }
 
   return distrib;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeAverageAmplitude(const Surface<Real>& surface) {
   return surface.mean();
 }
 /* -------------------------------------------------------------------------- */
 
 Complex SurfaceStatistics::computeAverageAmplitude(
     const SurfaceComplex<Real>& surface) {
   return surface.mean();
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeContactArea(const Surface<Real>& surface) {
   Real L = surface.getL();
   return SurfaceStatistics::computeContactAreaRatio(surface) * L * L;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SurfaceStatistics::computeContactAreaRatio(const Surface<Real>& surface) {
   UInt area = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for collapse(2), reduction(+ : area)
   for (UInt i = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       if (surface(i, j) > 1e-10)
         area++;
     }
   }
   return area / Real(n * n);
 }
 
 /* -------------------------------------------------------------------------- */
 
 Surface<Real>
 SurfaceStatistics::computeAutocorrelation(Surface<Real>& surface) {
   SurfaceComplex<Real> correlation_heights = computePowerSpectrum(surface);
   Surface<Real> acf(correlation_heights.size(), correlation_heights.getL());
   FFTransform<Real, 2>& plan =
       FFTPlanManager::get().createPlan(acf, correlation_heights);
   plan.backwardTransform();
   FFTPlanManager::get().destroyPlan(acf, correlation_heights);
   return acf;
 }
 /* -------------------------------------------------------------------------- */
 
 SurfaceComplex<Real>
 SurfaceStatistics::computePowerSpectrum(Surface<Real>& surface) {
   SurfaceComplex<Real> power_spectrum(surface.size(), surface.getL());
   FFTransform<Real, 2>& plan =
       FFTPlanManager::get().createPlan(surface, power_spectrum);
   plan.forwardTransform();
 
   UInt n = power_spectrum.dataSize();
   UInt factor = power_spectrum.size() * power_spectrum.size();
 
 #pragma omp parallel for
   for (UInt i = 0; i < n; ++i)
     power_spectrum(i) /= factor;
 
   power_spectrum.makeItRealBySquare();
   FFTPlanManager::get().destroyPlan(surface, power_spectrum);
   return power_spectrum;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 T SurfaceStatistics::computeAverage(const Surface<T>& surface, T epsilon) {
 
   Real avg = 0.0;
   UInt cpt = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+ : avg, cpt)
   for (UInt i = 0; i < n * n; i++) {
     Real val = surface(i);
     if (std::abs(val) > epsilon) {
       avg += val;
       ++cpt;
     }
   }
   avg /= cpt;
   return avg;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSkewness(const Surface<Real>& surface, Real avg,
                                         Real std) {
   Real skw = 0.0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+ : skw)
   for (UInt i = 0; i < n * n; i++) {
     Real tmp = surface(i) - avg;
     skw += (tmp * tmp * tmp);
   }
   skw /= (n * std * std * std);
   return skw;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeKurtosis(const Surface<Real>& surface, Real avg,
                                         Real std) {
 
   Real kur = 0.0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+ : kur)
   for (UInt i = 0; i < n * n; i++) {
     Real tmp = surface.size() - avg;
     kur += (tmp * tmp * tmp * tmp);
   }
   kur /= (n * std * std * std * std);
   return kur - 3;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeMaximum(const Surface<Real>& surface) {
   Real _max = -1. * std::numeric_limits<Real>::max();
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(max : _max)
   for (UInt i = 0; i < n * n; i++) {
     _max = std::max(surface(i), _max);
   }
   return _max;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeMinimum(const Surface<Real>& surface) {
   Real _min = std::numeric_limits<Real>::max();
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(min : _min)
   for (UInt i = 0; i < n * n; i++) {
     _min = std::min(surface(i), _min);
   }
   return _min;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeRMSSlope(const Surface<Real>& surface) {
 
   UInt n = surface.size();
   Surface<Real> grad(n, 1.);
 
   for (UInt i = 0, m = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       UInt k = i + n * j;
 
       //      xl = ((i==0)?(k+n-1):(k-1));
       UInt xh = ((i == (n - 1)) ? (k - n + 1) : (k + 1));
       //      yl = ((j==0)?(k+n*n-n):(k-n));
       UInt yh = ((j == (n - 1)) ? (k - n * n + n) : (k + n));
 
       Real gx = surface(xh) - surface(k);
       Real gy = surface(yh) - surface(k);
       grad(m++) = gx * gx + gy * gy;
     }
   }
 
   Real avg = SurfaceStatistics::computeAverage(grad);
   Real slp = sqrt(avg); /* According to paper Hyun PRE 70:026117 (2004) */
   Real scale_factor = surface.getL() / n;
   return slp / scale_factor;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSpectralRMSSlope(Surface<Real>& surface) {
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   UInt n = surface.size();
 
   Real rms_slopes = 0;
   // /!\ commented code was causing race conditions with dummy variable in
   // GridHermitian
 #pragma omp parallel for collapse(2), reduction(+ : rms_slopes)
   for (UInt i = 1; i < n / 2; ++i)
     for (UInt j = 1; j < n / 2; ++j) {
       rms_slopes += 2. * ((i * i) + (j * j)) * power(i, j).real();
       // rms_slopes += 1.*((i*i)+(j*j))*power(n-i,n-j).real();
       rms_slopes += 2. * ((i * i) + (j * j)) * power(n - i, j).real();
       // rms_slopes += 1.*((i*i)+(j*j))*power(i,n-j).real();
     }
 
       // external line (i or j == 0)
 #pragma omp parallel for reduction(+ : rms_slopes)
   for (UInt i = 1; i < n / 2; ++i) {
     rms_slopes += 1. * (i * i) * power(i, 0).real();
     rms_slopes += 1. * (i * i) * power(n - i, 0).real();
     rms_slopes += 2. * (i * i) * power(0, i).real();
     // rms_slopes += 1.*(i*i)*power(0,n-i).real();
   }
 
     // internal line (i or j == n/2)
 #pragma omp parallel for reduction(+ : rms_slopes)
   for (UInt i = 1; i < n / 2; ++i) {
     rms_slopes += 0.5 * (i * i + n * n / 4) * power(i, n / 2).real();
     rms_slopes += 0.5 * (i * i + n * n / 4) * power(n - i, n / 2).real();
     rms_slopes += 1. * (i * i + n * n / 4) * power(n / 2, i).real();
     // rms_slopes += 0.5*(i*i+n*n/4)*power(n/2,n-i).real();
   }
 
   {
     // i == n/2 j == n/2
     rms_slopes += 0.25 * (n * n / 2) * power(n / 2, n / 2).real();
     // i == 0 j == n/2
     rms_slopes += 1. * (n * n / 4) * power(0, n / 2).real();
     // i == n/2 j == 0
     rms_slopes += 1. * (n * n / 4) * power(n / 2, 0).real();
   }
   rms_slopes = 2. * M_PI / n * sqrt(rms_slopes);
   Real scale_factor = surface.getL() / n;
   return rms_slopes / scale_factor;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeStdev(const Grid<Real, 2>& surface,
                                      Real epsilon) {
   UInt n = surface.dataSize();
   UInt cpt = 0;
   Real std = 0.0;
   Real avg = 0.0;
 
 #pragma omp parallel for reduction(+ : avg, std, cpt)
   for (UInt i = 0; i < n; i++) {
     Real val = surface(i);
     if (val > epsilon) {
       avg += val;
       std += val * val;
       ++cpt;
     }
   }
   avg /= cpt;
   return sqrt(std / cpt - avg * avg);
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSpectralStdev(Surface<Real>& surface) {
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   Real rms_heights = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for collapse(2), reduction(+ : rms_heights)
   for (UInt i = 0; i < n; ++i)
     for (UInt j = 0; j < n; ++j) {
       rms_heights += power(i, j).real();
     }
   rms_heights = sqrt(rms_heights);
   return rms_heights;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSpectralMeanCurvature(Surface<Real>& surface) {
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   UInt n = surface.size();
 
   Real rms_curvatures = 0;
 #pragma omp parallel for collapse(2), reduction(+ : rms_curvatures)
   for (UInt i = 1; i < n / 2; ++i)
     for (UInt j = 1; j < n / 2; ++j) {
       Real coeff = ((i * i) + (j * j));
       coeff *= coeff;
       rms_curvatures += 1. * coeff * power(i, j).real();
       rms_curvatures += 1. * coeff * power(n - i, n - j).real();
       rms_curvatures += 1. * coeff * power(n - i, j).real();
       rms_curvatures += 1. * coeff * power(i, n - j).real();
     }
 
       // external line (i or j == 0)
 #pragma omp parallel for reduction(+ : rms_curvatures)
   for (UInt i = 1; i < n / 2; ++i) {
     Real coeff = i * i;
     coeff *= coeff;
     rms_curvatures += 4. / 3. * coeff * power(i, 0).real();
     rms_curvatures += 4. / 3. * coeff * power(n - i, 0).real();
     rms_curvatures += 4. / 3. * coeff * power(0, i).real();
     rms_curvatures += 4. / 3. * coeff * power(0, n - i).real();
   }
 
     // internal line (i or j == n/2)
 #pragma omp parallel for reduction(+ : rms_curvatures)
   for (UInt i = 1; i < n / 2; ++i) {
     Real coeff = (i * i + n * n / 4);
     coeff *= coeff;
     rms_curvatures += 0.5 * coeff * power(i, n / 2).real();
     rms_curvatures += 0.5 * coeff * power(n - i, n / 2).real();
     rms_curvatures += 0.5 * coeff * power(n / 2, i).real();
     rms_curvatures += 0.5 * coeff * power(n / 2, n - i).real();
   }
 
   {
     // i == n/2 j == n/2
     rms_curvatures +=
         0.25 * (n * n / 2) * (n * n / 2) * power(n / 2, n / 2).real();
     // i == 0 j == n/2
     rms_curvatures +=
         2. / 3. * (n * n / 4) * (n * n / 4) * power(0, n / 2).real();
     // i == n/2 j == 0
     rms_curvatures +=
         2. / 3. * (n * n / 4) * (n * n / 4) * power(n / 2, 0).real();
   }
 
   rms_curvatures *= 9. / 4.;
   rms_curvatures = M_PI / n / n * sqrt(rms_curvatures);
 
   return rms_curvatures;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computePerimeter(const Surface<Real>& surface,
                                          Real eps) {
   int perimeter = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for collapse(2), reduction(+ : perimeter)
   for (UInt i = 0; i < n; ++i) {
     for (UInt j = 0; j < n; ++j) {
       if (j > 0 && ((surface(i, (j - 1)) < eps && surface(i, j) > eps) ||
                     (surface(i, (j - 1)) > eps && surface(i, j) < eps)))
         perimeter++;
       if (i > 0 && ((surface(j, (i - 1)) < eps && surface(j, i) > eps) ||
                     (surface(j, (i - 1)) > eps && surface(j, i) < eps)))
         perimeter++;
     }
   }
 
   return perimeter / Real(n);
 }
 /* -------------------------------------------------------------------------- */
 Real SurfaceStatistics::computeSum(const Surface<Real>& surface) {
   Real sum = 0.;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+ : sum)
   for (UInt i = 0; i < n * n; ++i)
     sum += surface(i);
 
   return sum;
 }
 
 /* -------------------------------------------------------------------------- */
 std::vector<Real> SurfaceStatistics::computeMoments(Surface<Real>& surface) {
 
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   Real L = surface.getL();
   UInt n = surface.size();
   Real k_factor = 2. * M_PI / L;
 
   std::vector<std::vector<Real>> moment;
   moment.resize(6);
   for (int i = 0; i < 5; i++)
     moment[i].resize(5);
 
   for (UInt i = 0; i < n / 2; i++) {
     Real i0 = pow(Real(i) * k_factor, 0);
     Real i2 = pow(Real(i) * k_factor, 2);
     Real i4 = pow(Real(i) * k_factor, 4);
     Real im0 = pow(-Real(i) * k_factor, 0);
     Real im2 = pow(-Real(i) * k_factor, 2);
     Real im4 = pow(-Real(i) * k_factor, 4);
 
     for (UInt j = 0; j < n / 2; j++) {
       Real Phi = power(i, j).real();
 
       Real Phi_10;
       Real Phi_01;
       Real Phi_11;
       if (i == 0 && j == 0) {
         Phi_10 = 0;
         Phi_11 = 0;
         Phi_01 = 0;
       } else if (i == 0) {
         Phi_10 = 0;
         Phi_11 = 0;
         Phi_01 = power(i, (n - j)).real();
       } else if (j == 0) {
         Phi_10 = power(n - i, j).real();
         Phi_11 = 0;
         Phi_01 = 0;
       } else {
         Phi_10 = power(n - i, j).real();
         Phi_11 = power(n - i, (n - j)).real();
         Phi_01 = power(i, (n - j)).real();
       }
       Real j0 = pow(Real(j) * k_factor, 0);
       Real j2 = pow(Real(j) * k_factor, 2);
       Real j4 = pow(Real(j) * k_factor, 4);
       Real jm0 = pow(-Real(j) * k_factor, 0);
       Real jm2 = pow(-Real(j) * k_factor, 2);
       Real jm4 = pow(-Real(j) * k_factor, 4);
 
       moment[0][0] += i0 * j0 * Phi + im0 * j0 * Phi_10 + i0 * jm0 * Phi_01 +
                       im0 * jm0 * Phi_11;
       moment[0][2] += i0 * j2 * Phi + im0 * j2 * Phi_10 + i0 * jm2 * Phi_01 +
                       im0 * jm2 * Phi_11;
       moment[2][0] += i2 * j0 * Phi + im2 * j0 * Phi_10 + i2 * jm0 * Phi_01 +
                       im2 * jm0 * Phi_11;
       moment[2][2] += i2 * j2 * Phi + im2 * j2 * Phi_10 + i2 * jm2 * Phi_01 +
                       im2 * jm2 * Phi_11;
       moment[0][4] += i0 * j4 * Phi + im0 * j4 * Phi_10 + i0 * jm4 * Phi_01 +
                       im0 * jm4 * Phi_11;
       moment[4][0] += i4 * j0 * Phi + im4 * j0 * Phi_10 + i4 * jm0 * Phi_01 +
                       im4 * jm0 * Phi_11;
     }
   }
 
   std::vector<Real> mean_moments;
   // m2
   Real mean_m2 = 0.5 * (moment[0][2] + moment[2][0]);
   mean_moments.push_back(mean_m2);
   // m4
   Real mean_m4 = (3 * moment[2][2] + moment[0][4] + moment[4][0]) / 3.;
   mean_moments.push_back(mean_m4);
   return mean_moments;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeAnalyticFractalMoment(UInt order, UInt k1,
                                                      UInt k2, Real hurst,
                                                      Real A, Real L) {
   Real xi = k2 / k1;
   Real exp = -2. * hurst;
   Real _k1 = 2. * M_PI / L * k1;
 
   Real Cq;
   switch (order) {
   case 0:
     Cq = 2. * M_PI;
     break;
   case 2:
     Cq = M_PI;
     break;
   case 4:
     Cq = 3. / 4. * M_PI;
     break;
   default:
     SURFACE_FATAL("invalid order");
   }
   exp = order - 2. * hurst;
   Real res = A * Cq * pow(_k1, exp) / exp * (pow(xi, exp) - 1);
 
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
 
 #endif /* __SURFACE_STATISTICS_HH__ */
diff --git a/src/model/boussinesq.cpp b/src/model/boussinesq.cpp
index 75de5cc..3485511 100644
--- a/src/model/boussinesq.cpp
+++ b/src/model/boussinesq.cpp
@@ -1,163 +1,163 @@
 /**
  * @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 "boussinesq.hh"
 #include "influence.hh"
 #include "kelvin_helper.hh"
 #include "model.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 namespace detail {
 template <model_type type, typename T>
 class BoussinesqHelper {
   using boussinesq_t = T;
   using trait = model_type_traits<type>;
   static constexpr UInt dim = trait::dimension;
   static constexpr UInt bdim = trait::boundary_dimension;
 
 public:
   BoussinesqHelper(const T& boussinesq, Real y_3, Real cutoff)
       : boussinesq(boussinesq), 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,
             VectorProxy<Complex, dim>&& traction,
             VectorProxy<const Real, bdim>&& q) {
           // No Cutoff
           // if (-q.l2norm() * std::abs(y_3) < std::log(cutoff))
           //   return;
 
           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);
   }
 
   inline void applyConstantTerm(typename KelvinTrait<T>::out_t& out,
-                                typename KelvinTrait<T>::source_t& source,
-                                const influence::ElasticHelper<dim>& el) const {
+                                typename KelvinTrait<T>::source_t& /*source*/,
+                                const influence::ElasticHelper<dim>& /*el*/) const {
     out = 0;
   }
 
 protected:
   const T& boussinesq;
   const Real y_3, cutoff;
 };
 
 template <>
 inline void
 BoussinesqHelper<model_type::volume_2d, influence::Boussinesq<3, 1>>::
     applyConstantTerm(
         typename KelvinTrait<influence::Boussinesq<3, 1>>::out_t& out,
         typename KelvinTrait<influence::Boussinesq<3, 1>>::source_t& source,
         const influence::ElasticHelper<dim>& el) const {
   out = 0;
   out(0, 2) = -source(0) / el.mu;
   out(1, 2) = -source(1) / el.mu;
   out(2, 2) = -source(2) / (el.lambda + 2 * el.mu);
 }
 
 }  // namespace detail
 
 template <>
 Boussinesq<model_type::volume_2d, 0>::Boussinesq(Model* model) : parent(model) {
   this->initialize(trait::dimension, trait::dimension);
 }
 
 template <>
 Boussinesq<model_type::volume_2d, 1>::Boussinesq(Model* model) : parent(model) {
   this->initialize(trait::dimension, trait::dimension * trait::dimension);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Operator implementation */
 /* -------------------------------------------------------------------------- */
 template <model_type type, UInt derivative>
 void Boussinesq<type, derivative>::apply(GridBase<Real>& source,
                                          GridBase<Real>& out) const {
   Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus();
   influence::Boussinesq<trait::dimension, derivative> boussinesq(mu, nu);
 
   auto apply = [&](UInt i, decltype(this->source_buffers)& source_buffers,
                    decltype(this->disp_buffer)& out_buffer) {
     const Real L = this->model->getSystemSize().front();
     const UInt N = this->model->getDiscretization().front();
     const Real dl = L / (N - 1);
 
     out_buffer = 0;
 
     const Real xi = i * dl;
     detail::BoussinesqHelper<type, decltype(boussinesq)> helper(boussinesq, xi,
                                                                 1e-2);
     helper.applyIntegral(out_buffer, source_buffers.front(), this->wavevectors);
     typename detail::KelvinTrait<decltype(boussinesq)>::out_t out_fundamental(
         out_buffer(0));
     typename detail::KelvinTrait<decltype(boussinesq)>::source_t in_fundamental(
         source_buffers.front()(0));
 
     // Uniform shift in case of gradient computation
     helper.applyConstantTerm(out_fundamental, in_fundamental,
                              influence::ElasticHelper<trait::dimension>(mu, nu));
   };
 
   // the source is a Grid<Real, bdim> grid normally, so we gotta make sure we
   // can cast it into a Grid<Real, dim>
   GridView<Grid, Real, trait::boundary_dimension, trait::dimension> view(
       source, {}, -1);
   this->fourierApply(apply, view, out);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type, UInt derivative>
 void Boussinesq<type, derivative>::initialize(UInt source_components,
                                               UInt out_components) {
   // Copy horizontal sizes
   std::array<UInt, trait::boundary_dimension> sizes;
   std::copy(this->model->getDiscretization().begin() + 1,
             this->model->getDiscretization().end(), sizes.begin());
 
   auto hermitian_sizes =
       GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
           sizes);
   this->disp_buffer.setNbComponents(out_components);
   this->disp_buffer.resize(hermitian_sizes);
   this->source_buffers.resize(1);
   auto& source = this->source_buffers.front();
   source.setNbComponents(source_components);
   source.resize(hermitian_sizes);
 }
 
 __END_TAMAAS__
diff --git a/src/model/influence.hh b/src/model/influence.hh
index d7cef6f..87a67e6 100644
--- a/src/model/influence.hh
+++ b/src/model/influence.hh
@@ -1,498 +1,498 @@
 /**
  * @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"
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 namespace influence {
 
 namespace {
 const Complex I(0, 1);  /// < imaginary unit
 constexpr inline Real sign(bool upper) { return (upper) ? 1. : -1.; }
 }
 
 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,
+    VectorProxy<const Real, dim - 1>& /*q*/,
     StaticMatrix<Complex, ST, dim, dim>& f) const {
   Vector<Real, dim> e3{{{0, 0, 1}}};
   auto res = (f * e3);
   res *= -this->b;
   res(2) += 2. * f(2, 2);
   e3(2) *= -1 / (this->mu * this->b);  // << -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 shape_num, 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 shape_num, 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;
 // }
 
 namespace detail {
 template <UInt shape_num>
 struct LinearShape {
   using integrator = KelvinIntegrator<1>;
 
   template <bool domain, typename T>
   static inline T integrateG0(T& g0_mul, Real alpha, Real beta);
 
   template <bool domain, typename T>
   static inline T integrateG1(T& g1_mul, Real alpha, Real beta);
 };
 
 template <UInt shape_num>
 template <bool domain, typename T>
 inline T LinearShape<shape_num>::integrateG0(T& g0_mul, Real alpha,
                                              Real beta) {
   constexpr Int mul = (shape_num == 0) ? -1 : 1;
   const auto bounds = thrust::make_pair(std::min(-mul, 0), std::max(-mul, 0));
   return g0_mul * (integrator::g0<domain>(beta) *
                    (integrator::int_g0<domain>(bounds, alpha) +
                     mul * integrator::int_g0_z<domain>(bounds, alpha)));
 }
 
 template <UInt shape_num>
 template <bool domain, typename T>
 inline T LinearShape<shape_num>::integrateG1(T& g1_mul, Real alpha,
                                              Real beta) {
   constexpr Int mul = (shape_num == 0) ? -1 : 1;
   const auto bounds = thrust::make_pair(std::min(-mul, 0), std::max(-mul, 0));
   return g1_mul * (integrator::g1<domain>(beta) *
                        (integrator::int_g0<domain>(bounds, alpha) +
                         mul * integrator::int_g0_z<domain>(bounds, alpha)) +
                    integrator::g0<domain>(beta) *
                        (integrator::int_g1<domain>(bounds, alpha) +
                         mul * integrator::int_g1_z<domain>(bounds, alpha)));
 }
 }
 
 // Integrate against a linear shape function
 template <>
 template <Int yj_xi, UInt shape_num, 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;
 
   constexpr bool domain = (shape_num == 0) ? yj_xi >= 0 : yj_xi > 0;
 
   auto g0_mul = kelvin.template applyU0<domain>(q_, in);
   auto g1_mul = kelvin.template applyU1<domain>(q_, in);
 
   decltype(g0_mul) res;
 
   res = detail::LinearShape<shape_num>::template integrateG0<domain>(
       g0_mul, alpha, beta);
   res += detail::LinearShape<shape_num>::template integrateG1<domain>(
       g1_mul, alpha, beta);
 
   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/kelvin_helper.hh b/src/model/kelvin_helper.hh
index ca185cb..4b264c3 100644
--- a/src/model/kelvin_helper.hh
+++ b/src/model/kelvin_helper.hh
@@ -1,182 +1,182 @@
 /**
  * @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 KELVIN_HELPER_HH
 #define KELVIN_HELPER_HH
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "grid_hermitian.hh"
 #include "influence.hh"
 #include "model_type.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 namespace detail {
 
 /// Trait for kelvin local types
 template <typename T>
 struct KelvinTrait;
 
 template <UInt dim>
 struct KelvinTrait<influence::Kelvin<dim, 0>> {
   using source_t = VectorProxy<Complex, dim>;
   using out_t = VectorProxy<Complex, dim>;
 };
 
 template <UInt dim>
 struct KelvinTrait<influence::Kelvin<dim, 1>> {
   using source_t = MatrixProxy<Complex, dim, dim>;
   using out_t = VectorProxy<Complex, dim>;
 };
 
 template <UInt dim>
 struct KelvinTrait<influence::Kelvin<dim, 2>> {
   using source_t = MatrixProxy<Complex, dim, dim>;
   using out_t = MatrixProxy<Complex, dim, dim>;
 };
 
 template <UInt dim>
 struct KelvinTrait<influence::Boussinesq<dim, 0>> {
   using source_t = VectorProxy<Complex, dim>;
   using out_t = VectorProxy<Complex, dim>;
 };
 
 template <UInt dim>
 struct KelvinTrait<influence::Boussinesq<dim, 1>> {
   using source_t = VectorProxy<Complex, dim>;
   using out_t = MatrixProxy<Complex, dim, dim>;
 };
 
 /// Helper to apply integral representation on output
 template <model_type type, typename T>
 struct KelvinHelper {
   using kelvin_t = T;
   using trait = model_type_traits<type>;
   static constexpr UInt dim = trait::dimension;
   static constexpr UInt bdim = trait::boundary_dimension;
 
   KelvinHelper(const T& kelvin, Real cutoff) : kelvin(kelvin), cutoff(cutoff) {}
 
   template <Int yj_xi, UInt shape_num>
   inline void applyIntegral(GridHermitian<Real, bdim>& out,
                             GridHermitian<Real, bdim>& source, Real dij,
                             Real dl,
                             const Grid<Real, bdim>& wavevectors) const {
     Loop::stridedLoop(
         [&](typename KelvinTrait<T>::out_t&& out_local,
             typename KelvinTrait<T>::source_t&& source_local,
             VectorProxy<const Real, bdim>&& q) {
           // Cutoff
           if (-q.l2norm() * std::abs(dij) < std::log(cutoff))
             return;
           influence::KelvinIntegrator<1>::integrate<yj_xi, shape_num>(
               out_local, source_local, kelvin, q, dij, dl);
         },
         out, source, wavevectors);
   }
 
-  inline void applyFreeTerm(GridHermitian<Real, bdim>& out,
-                            GridHermitian<Real, bdim>& source,
-                            const Grid<Real, bdim>& wavevectors) const {}
+  inline void applyFreeTerm(GridHermitian<Real, bdim>& /*out*/,
+                            GridHermitian<Real, bdim>& /*source*/,
+                            const Grid<Real, bdim>& /*wavevectors*/) const {}
 
   const T& kelvin;
   const Real cutoff;
 };
 
 template <>
 inline void
 KelvinHelper<model_type::volume_2d, influence::Kelvin<3, 2>>::applyFreeTerm(
     GridHermitian<Real, bdim>& out, GridHermitian<Real, bdim>& source,
     const Grid<Real, bdim>& wavevectors) const {
   Loop::stridedLoop(
       [&](KelvinTrait<kelvin_t>::out_t&& out_local,
           KelvinTrait<kelvin_t>::source_t&& source_local,
           VectorProxy<const Real, bdim>&& q) {
         out_local += kelvin.applyDiscontinuityTerm(q, source_local);
       },
       out, source, wavevectors);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type, typename T>
 class LinearElement {
   using trait = model_type_traits<type>;
   static constexpr UInt bdim = trait::boundary_dimension;
 
 public:
   LinearElement(const thrust::pair<UInt, Real>& integration_node,
                 const T& kelvin,
                 std::vector<GridHermitian<Real, bdim>>& sources,
                 GridHermitian<Real, bdim>& out,
                 const Grid<Real, bdim>& wavevectors,
 		Real cutoff = 1e-2)
       : integration_node(integration_node), sources(sources), out(out),
         wavevectors(wavevectors), helper(kelvin, cutoff) {}
 
 private:
   template <UInt shape_num>
   void integrateNode(UInt element, Real yj, Real el_size) {
     static_assert(shape_num == 0 or shape_num == 1, "Shape number is invalid");
     const Real dij = yj - integration_node.second;
     const UInt node = element + shape_num;
 
     if (node > integration_node.first)
       helper.template applyIntegral<1, shape_num>(
           out, sources[node], dij, el_size, wavevectors);
     else if (node == integration_node.first) {
       helper.template applyIntegral<0, shape_num>(
           out, sources[node], dij, el_size, wavevectors);
 
       if (shape_num == 0) // apply free term only once
 	helper.applyFreeTerm(out, sources[node], wavevectors);
     }
 
     else
       helper.template applyIntegral<-1, shape_num>(
           out, sources[node], dij, el_size, wavevectors);
   }
 
 public:
   void integrateElement(UInt e, const thrust::pair<Real, Real>& positions) {
     const Real dl = std::abs(positions.second - positions.first);
     integrateNode<0>(e, positions.first, dl);
     integrateNode<1>(e, positions.second, dl);
   }
 
 private:
   thrust::pair<UInt, Real> integration_node;
   std::vector<GridHermitian<Real, bdim>>& sources;
   GridHermitian<Real, bdim>& out;
   const Grid<Real, bdim>& wavevectors;
   KelvinHelper<type, T> helper;
 };
 
 }  // namespace detail
 
 __END_TAMAAS__
 
 #endif  // KELVIN_HELPER_HH
diff --git a/src/model/mindlin.cpp b/src/model/mindlin.cpp
index 860ac71..3702fa7 100644
--- a/src/model/mindlin.cpp
+++ b/src/model/mindlin.cpp
@@ -1,185 +1,185 @@
 /**
  * @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"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 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>& surface_gradients,
                             GridHermitian<Real, bdim>& source_stresses,
                             const Grid<Real, bdim>& wavevectors) const {
     Loop::stridedLoop(
         [&](typename KelvinTrait<T>::out_t&& out_local,
             MatrixProxy<Complex, dim, dim>&& surface_gradient,
             MatrixProxy<Complex, dim, dim>&& source_stress,
             VectorProxy<const Real, bdim>&& q) {
           // No cutoff
           // if (-q.l2norm() * std::abs(y_3) < std::log(cutoff))
           //   return;
 
           Vector<Real, dim> normal{{{0, 0, 1}}};
           auto traction =
               el(surface_gradient) * normal;   // computing strains from gradu
           traction -= source_stress * normal;  // subtracting source stresses
           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, surface_gradients, source_stresses, wavevectors);
   }
 
   inline void applyConstantTerm(typename KelvinTrait<T>::out_t& out,
-                                Vector<Complex, dim>& source,
-                                const influence::ElasticHelper<dim>& el) const {
+                                Vector<Complex, dim>& /*source*/,
+                                const influence::ElasticHelper<dim>& /*el*/) const {
     out = 0;
   }
 
 protected:
   const T& boussinesq;
   const influence::ElasticHelper<dim>& el;
   const Real y_3, cutoff;
 };
 
 template <>
 inline void
 MindlinBoussinesqHelper<model_type::volume_2d, influence::Boussinesq<3, 1>>::
     applyConstantTerm(
         typename KelvinTrait<influence::Boussinesq<3, 1>>::out_t& out,
         Vector<Complex, dim>& source,
         const influence::ElasticHelper<dim>& el) const {
   out = 0;
   out(2, 0) = -source(0) / el.mu;
   out(2, 1) = -source(1) / el.mu;
   out(2, 2) = -source(2) / (el.lambda + 2 * el.mu);
 }
 }  // namespace detail
 
 template <model_type type, UInt order>
 void Mindlin<type, order>::applyIf(GridBase<Real>& source, GridBase<Real>& out,
                                    filter_t pred) 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 node, decltype(this->source_buffers)& source_buffers,
                    decltype(this->disp_buffer)& out_buffer) {
     const Real L = this->model->getSystemSize().front();
     const UInt N = this->model->getDiscretization().front();
     const Real dl = L / (N - 1);
     const Real xi = node * dl;
 
     // this should only be used when node == xi == 0
     detail::LinearElement<type, decltype(kelvin_strain)> strain_element{
         thrust::make_pair(node, xi), kelvin_strain, source_buffers,
         surface_strains, this->wavevectors};
     detail::LinearElement<type, decltype(kelvin)> element{
         thrust::make_pair(node, xi), kelvin, source_buffers, out_buffer,
         this->wavevectors};
 
     out_buffer = 0;
 
     // Simple condition to ensure correct integration
     auto filter_out_elem = [&pred](UInt e) {
       return not pred(e) and not pred(e + 1);
     };
 
     /// Compute surface strains only once (dirty)
     if (node == 0) {
       surface_strains = 0;
       for (UInt e : Loop::range(N - 1)) {
         if (filter_out_elem(e))
           continue;
 
         const auto node_positions = thrust::make_pair(e * dl, (e + 1) * dl);
         strain_element.integrateElement(e, node_positions);
       }
     }
 
     // Apply Kelvin operator
     for (UInt e : Loop::range(N - 1)) {
       if (filter_out_elem(e))
         continue;
 
       const auto node_positions = thrust::make_pair(e * dl, (e + 1) * dl);
       element.integrateElement(e, node_positions);
     }
 
     // Correcting for the tractions on the surface
     detail::MindlinBoussinesqHelper<type, decltype(boussinesq)> helper(
         boussinesq, elasticity, xi, 1e-3);
     helper.applyIntegral(out_buffer, surface_strains, source_buffers.front(),
                          this->wavevectors);
 
     typename detail::KelvinTrait<decltype(boussinesq)>::out_t out_fundamental(
         out_buffer(0));
     typename detail::KelvinTrait<decltype(kelvin)>::source_t source_fundamental(
         source_buffers[node](0));
 
     Vector<Real, trait::dimension> e3;
     e3 = 0;
     e3(trait::dimension - 1) = -1;
     auto in_fundamental = source_fundamental * e3;
 
     // Uniform shift in case of gradient computation
     helper.applyConstantTerm(out_fundamental, in_fundamental, elasticity);
   };
 
   this->fourierApplyIf(apply, source, out, pred);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation                                                     */
 /* -------------------------------------------------------------------------- */
 template class Mindlin<model_type::volume_2d, 3>;
 template class Mindlin<model_type::volume_2d, 4>;
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
diff --git a/src/solvers/contact_solver.hh b/src/solvers/contact_solver.hh
index ec9e3be..219329c 100644
--- a/src/solvers/contact_solver.hh
+++ b/src/solvers/contact_solver.hh
@@ -1,77 +1,77 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016-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 __CONTACT_SOLVER_HH__
 #define __CONTACT_SOLVER_HH__
 /* -------------------------------------------------------------------------- */
 #include "model.hh"
 #include "tamaas.hh"
 #include "meta_functional.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 
 class ContactSolver {
 public:
   /// Constructor
   ContactSolver(Model& model, const GridBase<Real>& surface, Real tolerance);
   /// Destructor
   virtual ~ContactSolver();
 
 public:
   /// Print state of solve
   virtual void printState(UInt iter, Real cost_f, Real error);
   /// Set maximum number of iterations
   void setMaxIterations(UInt n) { max_iterations = n; }
   /// Set dump_frequency
   void setDumpFrequency(UInt n) { dump_frequency = n; }
   /// Add term to functional
   void addFunctionalTerm(functional::Functional* func);
 
   // Virtual functions
 public:
-  virtual Real solve(Real load) {
+  virtual Real solve(Real /*load*/) {
     TAMAAS_EXCEPTION("solve() not implemented");
   }
 
 public:
   GridBase<Real>& getSurface() { return surface; }
   Model& getModel() { return model; }
 
 protected:
   Model& model;
   GridBase<Real> surface;
   std::unique_ptr<GridBase<Real>> _gap = nullptr;
   functional::MetaFunctional functional;
 
   Real tolerance;
   UInt max_iterations = 1000;
   Real surface_stddev;
   UInt dump_frequency = 100;
 };
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 #endif  // __CONTACT_SOLVER_HH__
diff --git a/src/solvers/ep_solver.cpp b/src/solvers/ep_solver.cpp
index be84121..40b3938 100644
--- a/src/solvers/ep_solver.cpp
+++ b/src/solvers/ep_solver.cpp
@@ -1,46 +1,46 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016-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 "model_type.hh"
 #include "ep_solver.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 EPSolver::EPSolver(Residual& residual)
     : _residual(residual) {
   auto& model = residual.getModel();
   auto dim = model.getDiscretization().size();
   _x = allocateGrid<false, Real>(model.getType(), model.getDiscretization(),
                                  dim * dim);
 }
 
 /* -------------------------------------------------------------------------- */
 void EPSolver::updateState() {
-  _residual.updateState(*_x);
+  _residual.updateState(getStrainIncrement());
 }
 
 }  // namespace tamaas