diff --git a/src/model/model_factory.cpp b/src/model/model_factory.cpp
index 66f34cf..fa939cd 100644
--- a/src/model/model_factory.cpp
+++ b/src/model/model_factory.cpp
@@ -1,145 +1,145 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "model_factory.hh"
+#include "dcfft.hh"
 #include "materials/isotropic_hardening.hh"
 #include "model_template.hh"
-#include "dcfft.hh"
 #include <functional>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 template <typename Abstract, template <model_type> class Concrete,
           class... Args>
 decltype(auto) createFromModelType(model_type type, Args&&... args) {
   using namespace std::placeholders;
   return model_type_dispatch(
       std::bind(
           [](auto&& _, Args&&... args) -> std::unique_ptr<Abstract> {
             constexpr auto type = std::decay_t<decltype(_)>::value;
             return std::make_unique<Concrete<type>>(
                 std::forward<Args>(args)...);
           },
           _1, std::forward<Args>(args)...),
       type);
 }
 
 std::unique_ptr<Model>
 ModelFactory::createModel(model_type type, const std::vector<Real>& system_size,
                           const std::vector<UInt>& discretization) {
   auto model = createFromModelType<Model, ModelTemplate>(type, system_size,
                                                          discretization);
 
   const auto& shape = model->getDiscretization();
   if (std::any_of(std::begin(shape), std::end(shape),
                   [](auto&& v) { return v == 0; }))
     TAMAAS_EXCEPTION("FFT data partition gave no data on ["
                      << mpi::rank() << "|" << mpi::size() << "]");
 
   return model;
 }
 
 std::unique_ptr<Model> ModelFactory::createModel(const Model& model) {
   auto copy = createModel(model.getType(), model.getSystemSize(),
                           model.getGlobalDiscretization());
 
   copy->setElasticity(model.getYoungModulus(), model.getPoissonRatio());
 
   // Register and copy all fields
   for (auto&& pair : model.fields_map()) {
     boost::apply_visitor(
         [&copy, &pair](auto&& grid_ptr) {
           auto&& grid = *grid_ptr;
           dimension_dispatch(
               [&copy, &grid, &pair](auto&& _) {
                 using T = typename std::decay_t<decltype(grid)>::value_type;
                 constexpr UInt dim = std::decay_t<decltype(_)>::value;
                 const auto* ptr = dynamic_cast<const Grid<T, dim>*>(&grid);
                 if (ptr != nullptr) {
                   auto new_f = std::make_unique<Grid<T, dim>>(
                       ptr->sizes(), ptr->getNbComponents());
                   *new_f = *ptr;
                   (*copy)[pair.first] = std::move(new_f);
                 }
               },
               grid.getDimension());
         },
         pair.second);
   }
 
   // TODO: Registrer and copy integral operators
 
   return copy;
 }
 
 std::unique_ptr<Residual>
 ModelFactory::createResidual(Model& model, Real sigma_y, Real hardening) {
   if (model.getType() != model_type::volume_2d)
     TAMAAS_EXCEPTION("Cannot instanciate model: " << model);
 
   // Make elastic-plastic material
   auto mat = std::make_unique<IsotropicHardening>(&model, sigma_y, hardening);
   return std::make_unique<Residual>(model, std::move(mat));
 }
 
 void ModelFactory::registerVolumeOperators(Model& m) {
   if (m.getType() != model_type::volume_2d)
     TAMAAS_EXCEPTION("Registering volume operators not supported on "
                      << m.getType());
 
   constexpr auto type = model_type::volume_2d;
   m.registerIntegralOperator<Mindlin<type, 2>>("mindlin_gradient");
   m.registerIntegralOperator<Boussinesq<type, 1>>("boussinesq_gradient");
   m.registerIntegralOperator<Mindlin<type, 1>>("mindlin");
   m.registerIntegralOperator<Boussinesq<type, 0>>("boussinesq");
 }
 
 void ModelFactory::registerNonPeriodic(Model& m, std::string name) {
   if (m.getType() != model_type::basic_2d)
     TAMAAS_EXCEPTION("Registering non-periodic operator not supported on "
                      << m.getType());
 
   m.registerIntegralOperator<DCFFT>(name);
 }
 
 void ModelFactory::setIntegrationMethod(IntegralOperator& op,
                                         integration_method method,
                                         Real cutoff) {
 #define CAST(derivative)                                                       \
   do {                                                                         \
     auto* casted =                                                             \
         dynamic_cast<Kelvin<model_type::volume_2d, derivative>*>(&op);         \
     if (casted != nullptr) {                                                   \
       casted->setIntegrationMethod(method, cutoff);                            \
       return;                                                                  \
     }                                                                          \
   } while (0)
 
   CAST(0);
   CAST(1);
   CAST(2);
 
 #undef CAST
 }
 
 }  // namespace tamaas