diff --git a/src/model/model.cpp b/src/model/model.cpp
index 65f6ccc..dad6c59 100644
--- a/src/model/model.cpp
+++ b/src/model/model.cpp
@@ -1,191 +1,188 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "model.hh"
 #include "be_engine.hh"
 #include "logger.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 Model::Model(std::vector<Real> system_size, std::vector<UInt> discretization)
     : system_size(std::move(system_size)),
       discretization(std::move(discretization)) {}
 
 /* -------------------------------------------------------------------------- */
 
 Model::~Model() = default;
 
 /* -------------------------------------------------------------------------- */
 
 void Model::setElasticity(Real E_, Real nu_) {
   setYoungModulus(E_);
   setPoissonRatio(nu_);
   updateOperators();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::applyElasticity(GridBase<Real>& stress,
                             const GridBase<Real>& strain) const {
   operators.at("hooke")->apply(const_cast<GridBase<Real>&>(strain), stress);
 }
 
 /* -------------------------------------------------------------------------- */
 
-Real Model::getHertzModulus() const { return E / (1 - nu * nu); }
-
-/* -------------------------------------------------------------------------- */
-
 GridBase<Real>& Model::getTraction() { return getField("traction"); }
-
 const GridBase<Real>& Model::getTraction() const {
   return getField("traction");
 }
 
 /* -------------------------------------------------------------------------- */
 
 GridBase<Real>& Model::getDisplacement() { return getField("displacement"); }
-
 const GridBase<Real>& Model::getDisplacement() const {
   return getField("displacement");
 }
 
 /* -------------------------------------------------------------------------- */
 
 const std::vector<Real>& Model::getSystemSize() const { return system_size; }
 
 /* -------------------------------------------------------------------------- */
 
 const std::vector<UInt>& Model::getDiscretization() const {
   return discretization;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void Model::solveNeumann() {
   engine->registerNeumann();
   engine->solveNeumann(getField("traction"), getField("displacement"));
 }
 
 void Model::solveDirichlet() {
   engine->registerDirichlet();
   engine->solveDirichlet(getField("displacement"), getField("traction"));
 }
 
 /* -------------------------------------------------------------------------- */
 IntegralOperator* Model::getIntegralOperator(const std::string& name) {
   return operators.at(name).get();
 }
 
 /* -------------------------------------------------------------------------- */
 std::vector<std::string> Model::getIntegralOperators() const {
   std::vector<std::string> keys;
   keys.reserve(operators.size());
 
   std::transform(operators.begin(), operators.end(), std::back_inserter(keys),
                  [](auto&& pair) { return pair.first; });
   return keys;
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::updateOperators() {
   for (auto& op : operators)
     op.second->updateFromModel();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumper(std::shared_ptr<ModelDumper> dumper) {
   this->dumpers.push_back(std::move(dumper));
 }
 
 void Model::dump() const {
   for (auto& dumper : dumpers)
     if (dumper)
       dumper->dump(*this);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::registerField(const std::string& name,
                           std::shared_ptr<GridBase<Real>> field) {
   fields[name] = std::move(field);
 }
 
 const GridBase<Real>& Model::getField(const std::string& name) const try {
   return *fields.at(name);
 } catch (std::out_of_range& e) {
   Logger().get(LogLevel::warning)
       << "Field " << name << " not registered in model\n";
   throw e;
 }
 
 GridBase<Real>& Model::getField(const std::string& name) try {
   return *fields.at(name);
 } catch (std::out_of_range& e) {
   Logger().get(LogLevel::warning)
       << "Field " << name << " not registered in model\n";
   throw e;
 }
 
 std::vector<std::string> Model::getFields() const {
   std::vector<std::string> keys;
   keys.reserve(fields.size());
 
   std::transform(fields.begin(), fields.end(), std::back_inserter(keys),
                  [](auto&& pair) { return pair.first; });
   return keys;
 }
 
 /* -------------------------------------------------------------------------- */
 std::ostream& operator<<(std::ostream& o, const Model& _this) {
   o << "Model<" << _this.getType() << "> (E = " << _this.getYoungModulus()
     << ", nu = " << _this.getPoissonRatio() << ")\n";
 
   auto out_collec = [&o](auto&& collec) {
     std::for_each(collec.begin(), collec.end() - 1,
                   [&o](const auto& x) { o << x << ", "; });
     o << collec.back();
   };
 
   // Printing domain size
   o << "  - domain = [";
   out_collec(_this.getSystemSize());
   o << "]\n";
 
   // Printing discretization
   o << "  - discretization = [";
   out_collec(_this.getDiscretization());
   o << "]\n";
 
   // Print fields
   o << "  - registered fields = [";
   out_collec(_this.getFields());
   o << "]\n";
 
   o << "  - registered operators = [";
   out_collec(_this.getIntegralOperators());
-  o << "]\n";
+  o << "]";
+
+  if (_this.dumpers.size())
+    o << "\n  - " << _this.dumpers.size() << " registered dumpers";
 
   return o;
 }
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
diff --git a/src/model/model.hh b/src/model/model.hh
index 36930ad..cb12a3f 100644
--- a/src/model/model.hh
+++ b/src/model/model.hh
@@ -1,195 +1,196 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef MODEL_HH
 #define MODEL_HH
 /* -------------------------------------------------------------------------- */
 #include "be_engine.hh"
 #include "grid_base.hh"
 #include "integral_operator.hh"
 #include "model_dumper.hh"
 #include "model_type.hh"
 #include "tamaas.hh"
 #include <algorithm>
 #include <memory>
 #include <unordered_map>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * @brief Model containing pressure and displacement
  * This class is a container for the model fields. It is supposed to be
  * dimension agnostic, hence the GridBase members.
  */
 class Model {
 protected:
   /// Constructor
   Model(std::vector<Real> system_size, std::vector<UInt> discretization);
 
 public:
   /// Destructor
   virtual ~Model();
 
 public:
   /// Set elasticity parameters
   void setElasticity(Real E, Real nu);
   /// Get Hertz contact modulus
-  Real getHertzModulus() const;
+  Real getHertzModulus() const { return E / (1 - nu * nu); }
   /// Get Young's modulus
   Real getYoungModulus() const { return E; }
   /// Get Poisson's ratio
   Real getPoissonRatio() const { return nu; }
   /// Get shear modulus
   Real getShearModulus() const { return E / (2 * (1 + nu)); }
   /// Set Young's modulus
   void setYoungModulus(Real E_) {
     if (E_ < 0)
       TAMAAS_EXCEPTION("Elastic modulus should be positive");
     this->E = E_;
     updateOperators();
   }
 
   /// Set Poisson's ratio
   void setPoissonRatio(Real nu_) {
     if (nu_ > 0.5 or nu_ <= -1)
       TAMAAS_EXCEPTION("Poisson ratio should be in ]-1, 0.5]");
     this->nu = nu_;
     updateOperators();
   }
 
-public:
-  void applyElasticity(GridBase<Real>& stress,
-                       const GridBase<Real>& strain) const;
-
 public:
   /// Get pressure
   GridBase<Real>& getTraction();
   /// Get pressure
   const GridBase<Real>& getTraction() const;
   /// Get displacement
   GridBase<Real>& getDisplacement();
   /// Get displacement
   const GridBase<Real>& getDisplacement() const;
   /// Get model type
   virtual model_type getType() const = 0;
   /// Get system physical size
   const std::vector<Real>& getSystemSize() const;
   /// Get boundary system physical size
   virtual std::vector<Real> getBoundarySystemSize() const = 0;
   /// Get discretization
   const std::vector<UInt>& getDiscretization() const;
   /// Get discretization of global MPI system
   virtual std::vector<UInt> getGlobalDiscretization() const = 0;
   /// Get boundary discretization
   virtual std::vector<UInt> getBoundaryDiscretization() const = 0;
   /// Get boundary element engine
   BEEngine& getBEEngine() {
     TAMAAS_ASSERT(engine, "BEEngine was not initialized");
     return *engine;
   }
 
+public:
+  /// Apply Hooke's law
+  void applyElasticity(GridBase<Real>& stress,
+                       const GridBase<Real>& strain) const;
+
   /// Solve Neumann problem using default neumann operator
   void solveNeumann();
   /// Solve Dirichlet problem using default dirichlet operator
   void solveDirichlet();
 
 public:
   /// Register a new integral operator
   template <typename Operator>
   IntegralOperator* registerIntegralOperator(const std::string& name) {
     Logger().get(LogLevel::debug)
         << TAMAAS_DEBUG_MSG("registering operator " + name);
     operators[name] = std::make_unique<Operator>(this);
     return operators[name].get();
   }
 
   /// Get a registerd integral operator
   IntegralOperator* getIntegralOperator(const std::string& name);
   /// Get list of integral operators
   std::vector<std::string> getIntegralOperators() const;
   /// Get operators mapcar
   const auto& getIntegralOperatorsMap() const { return operators; }
   /// Tell operators to update their cache
   void updateOperators();
 
 public:
   /// Register a field
   void registerField(const std::string& name,
                      std::shared_ptr<GridBase<Real>> field);
   /// Get a field
   const GridBase<Real>& getField(const std::string& name) const;
   /// Get a non-const field
   GridBase<Real>& getField(const std::string& name);
-
   /// Get fields
   std::vector<std::string> getFields() const;
-
   /// Get fields map
   const auto& getFieldsMap() const { return fields; }
 
 public:
   /// Set the dumper object
   void addDumper(std::shared_ptr<ModelDumper> dumper);
   /// Dump the model
   void dump() const;
 
+  friend std::ostream& operator<<(std::ostream& o, const Model& _this);
+
 protected:
   Real E = 1, nu = 0;
   std::vector<Real> system_size;
   std::vector<UInt> discretization;
   std::unique_ptr<BEEngine> engine = nullptr;
   std::unordered_map<std::string, std::shared_ptr<IntegralOperator>> operators;
   std::unordered_map<std::string, std::shared_ptr<GridBase<Real>>> fields;
   std::vector<std::shared_ptr<ModelDumper>> dumpers;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Output model to stream */
 /* -------------------------------------------------------------------------- */
 std::ostream& operator<<(std::ostream& o, const Model& _this);
 
 /* -------------------------------------------------------------------------- */
 /* Simpler grid allocation */
 /* -------------------------------------------------------------------------- */
 template <bool boundary, typename T>
 std::unique_ptr<GridBase<T>> allocateGrid(Model& model) {
   return allocateGrid<boundary, T>(
       model.getType(), (boundary) ? model.getBoundaryDiscretization()
                                   : model.getDiscretization());
 }
 
 template <bool boundary, typename T>
 std::unique_ptr<GridBase<T>> allocateGrid(Model& model, UInt nc) {
   return allocateGrid<boundary, T>(model.getType(),
                                    (boundary)
                                        ? model.getBoundaryDiscretization()
                                        : model.getDiscretization(),
                                    nc);
 }
 
 }  // namespace tamaas
 
 #endif  // MODEL_HH