diff --git a/CHANGELOG.md b/CHANGELOG.md
index f9d2f30..725a51b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,186 +1,187 @@
 # Changelog
 All notable changes to this project will be documented in this file.
 
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## Unreleased
 
 ### Added
 
 - Added `read()` method to dumpers to create a model from a dump file
 - `getClusters()` can be called in MPI contact with partial contact maps
 - Added a JSON encoder class for models and a JSON dumper
 
 ### Fixed
 
 - Fixed a host of dump read/write issues when model type was not `volume_*d`.
   Dumper tests are now streamlined and systematic.
+- Fixed a bug where `Model::solveDirichlet` would not compute correctly
 
 ## v2.2.2 -- 2021-04-02
 
 ### Added
 
 - Entry-point `tamaas` defines a grouped CLI for `examples/pipe_tools`. Try
   executing `tamaas surface -h` from the command-line!
 
 ### Changed
 
 - `CXXFLAGS` are now passed to the linker
 - Added this changelog
 - Using absolute paths for environmental variables when running `scons test`
 - Reorganized documentation layout
 - Gave the build system a facelift (docs are now generated directly with SCons
   instead of a Makefile)
 
 ### Deprecated
 
 - Python 2 support is discontinued. Version `v2.2.1` is the last PyPi build with
   a Python 2 wheel.
 - The scripts in `examples/pipe_tools` have been replaced by the `tamaas` command
 
 ### Fixed
 
 - `UVWDumper` no longer imports `mpi4py` in sequential
 - Compiling with different Thrust/FFTW backends
 
 
 ## v2.2.1 -- 2021-03-02
 
 ### Added
 
 - Output registered fields and dumpers in `print(model)`
 - Added `operator[]` to the C++ model class (for fields)
 - Added `traction` and `displacement` properties to Python model bindings
 - Added `operators` property to Python model bindings, which provides a
   dict-like access to registered operators
 - Added `shape` and `spectrum` to properties to Python surface generator
   bindings
 - Surface generator constructor accepts surface global shape as argument
 - Choice of FFTW thread model
 
 ### Changed
 
 - Tests use `/tmp` for temporary files
 - Updated dependency versions (Thrust, Pybind11)
 
 ### Deprecated
 
 - Most `get___()` and `set___()` in Python bindings have been deprecated. They
   will generate a `DeprecationWarning`.
 
 ### Removed
 
 - All legacy code
 
 
 ## v2.2.0 -- 2020-12-31
 
 ### Added
 
 - More accurate function for computation of contact area
 - Function to compute deviatoric of tensor fields
 - MPI implementation
 - Convenience `hdf5toVTK` function
 - Readonly properties `shape`, `global_shape`, `boundary_shape` on model to give
   shape information
 
 ### Changed
 
 - Preprocessor defined macros are prefixed with `TAMAAS_`
 - Moved `tamaas.to_voigt` to `tamaas.compute.to_voigt`
 
 ### Fixed
 
 - Warning about deprecated constructors with recent GCC versions
 - Wrong computation of grid strides
 - Wrong computation of grid sizes in views
 
 
 ## v2.1.4 -- 2020-08-07
 
 ### Added
 
 - Possibility to generate a static `libTamaas`
 - C++ implementation of DFSANE solver
 - Allowing compilation without OpenMP
 
 ### Changed
 
 - NetCDF dumper writes frames to a single file
 
 ### Fixed
 
 - Compatibility with SCons+Python 3
 
 ## v2.1.3 -- 2020-07-27
 
 ### Added
 
 - Version number to `TamaasInfo`
 
 ### Changed
 
 - Prepending root directory when generating archive
 
 
 ## v2.1.2 -- 2020-07-24
 
 This release changes some core internals related to discrete Fourier transforms
 for future MPI support.
 
 ### Added
 
 - Caching `CXXFLAGS` in SCons build
 - SCons shortcut to create code archive
 - Test of the elastic-plastic contact solver
 - Paraview data dumper (`.pvd` files)
 - Compression for UVW dumper
 - `__contains__` and `__iter__` Python bindings of model
 - Warning message of possible overflow in Kelvin
 
 ### Changed
 
 - Simplified `tamaas_info.cpp`, particularly the diff part
 - Using a new class `FFTEngine` to manage discrete Fourier transforms. Plans are
   re-used as much as possible with different data with the same shape. This is
   in view of future MPI developments
 - Redirecting I/O streams in solve functions so they can be used from Python
   (e.g. in Jupyter notebooks)
 - Calling `initialize()` and `finalize()` is no longer necessary
 
 ### Fixed
 
 - Convergence issue with non-linear solvers
 - Memory error in volume potentials
 
 
 ## v2.1.1 -- 2020-04-22
 
 ### Added
 
 - SCons shortcut to run tests
 
 ### Fixed
 
 - Correct `RPATH` for shared libraries
 - Issues with SCons commands introduced in v2.1.0
 - Tests with Python 2.7
 
 
 ## v2.1.0 -- 2020-04-17
 
 ### Added
 
 - SCons shortcuts to build/install Tamaas and its components
 - Selection of integration method for Kelvin operator
 - Compilation option to remove the legacy part of Tamaas
 - NetCDF dumper
 
 ### Fixed
 
 - Link bug with clang
 - NaNs in Kato saturated solver
 
 ## v2.0.0 -- 2019-11-11
 
 First public release. Contains relatively mature elastic-plastic contact code.
diff --git a/src/model/model.cpp b/src/model/model.cpp
index 16f2970..b358cf2 100644
--- a/src/model/model.cpp
+++ b/src/model/model.cpp
@@ -1,184 +1,184 @@
 /**
  *  @file
  *  LICENSE
  *
  *  Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "model.hh"
 #include "be_engine.hh"
 #include "logger.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 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);
 }
 
 /* -------------------------------------------------------------------------- */
 GridBase<Real>& Model::getTraction() { return getField("traction"); }
 const GridBase<Real>& Model::getTraction() const { return (*this)["traction"]; }
 
 /* -------------------------------------------------------------------------- */
 GridBase<Real>& Model::getDisplacement() { return getField("displacement"); }
 const GridBase<Real>& Model::getDisplacement() const {
   return (*this)["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(getTraction(), getDisplacement());
 }
 
 void Model::solveDirichlet() {
   engine->registerDirichlet();
-  engine->solveDirichlet(getTraction(), getDisplacement());
+  engine->solveDirichlet(getDisplacement(), getTraction());
 }
 
 /* -------------------------------------------------------------------------- */
 IntegralOperator* Model::getIntegralOperator(const std::string& name) const {
   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;
 }
 
 const GridBase<Real>& Model::operator[](const std::string& name) const {
   return getField(name);
 }
 
 GridBase<Real>& Model::operator[](const std::string& name) {
   return getField(name);
 }
 
 /* -------------------------------------------------------------------------- */
 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";
 
   if (mpi::size() > 1) {
     o << "  - global discretization = [";
     out_collec(_this.getGlobalDiscretization());
     o << "]\n";
   }
 
   // Print fields
   o << "  - registered fields = [";
   out_collec(_this.getFields());
   o << "]\n";
 
   o << "  - registered operators = [";
   out_collec(_this.getIntegralOperators());
   o << "]";
 
   if (_this.dumpers.size())
     o << "\n  - " << _this.dumpers.size() << " registered dumpers";
 
   return o;
 }
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas