diff --git a/src/common/ccoord_operations.hh b/src/common/ccoord_operations.hh
index 28f9d51..4e502aa 100644
--- a/src/common/ccoord_operations.hh
+++ b/src/common/ccoord_operations.hh
@@ -1,320 +1,320 @@
 /**
  * @file   ccoord_operations.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   29 Sep 2017
  *
  * @brief  common operations on pixel addressing 
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre 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
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include <functional>
 #include <numeric>
 #include <vector>
 #include <utility>
 
 #include <Eigen/Dense>
 
 #include "common/common.hh"
 #include "common/iterators.hh"
 
 #ifndef CCOORD_OPERATIONS_H
 #define CCOORD_OPERATIONS_H
 
 namespace muSpectre {
 
   namespace CcoordOps {
 
     namespace internal {
       //! simple helper returning the first argument and ignoring the second
       template <typename T>
       constexpr T ret(T val, size_t /*dummy*/) {return val;}
 
       //! helper to build cubes
       template <Dim_t Dim, typename T, size_t... I>
       constexpr std::array<T, Dim> cube_fun(T val, std::index_sequence<I...>) {
         return std::array<T, Dim>{ret(val, I)...};
       }
 
       //! computes hermitian size according to FFTW
       template <Dim_t Dim, size_t... I>
       constexpr Ccoord_t<Dim> herm(const Ccoord_t<Dim> & full_sizes,
                                    std::index_sequence<I...>) {
         return Ccoord_t<Dim>{full_sizes[I]..., full_sizes.back()/2+1};
       }
 
       //! compute the stride in a direction of a row-major grid
       template <Dim_t Dim>
       constexpr Dim_t stride(const Ccoord_t<Dim> & sizes,
                              const size_t index) {
         static_assert(Dim > 0, "only for positive numbers of dimensions");
 
         auto const diff{Dim - 1 - Dim_t(index)};
         Dim_t ret_val{1};
         for (Dim_t i{0}; i < diff; ++i) {
           ret_val *= sizes[Dim-1-i];
         }
         return ret_val;
       }
 
       //! get all strides from a row-major grid (helper function)
       template <Dim_t Dim, size_t... I>
       constexpr Ccoord_t<Dim> compute_strides(const Ccoord_t<Dim> & sizes,
                                               std::index_sequence<I...>) {
         return Ccoord_t<Dim>{stride<Dim>(sizes, I)...};
       }
     }  // internal
 
     //----------------------------------------------------------------------------//
     //! returns a grid of equal resolutions in each direction
     template <size_t dim, typename T>
     constexpr std::array<T, dim> get_cube(T size) {
       return internal::cube_fun<dim>(size, std::make_index_sequence<dim>{});
     }
 
     /* ---------------------------------------------------------------------- */
     //! returns the hermition grid to correcsponding to a full grid
     template <size_t dim>
     constexpr Ccoord_t<dim> get_hermitian_sizes(Ccoord_t<dim> full_sizes) {
       return internal::herm<dim>(full_sizes, std::make_index_sequence<dim-1>{});
     }
 
     /* ---------------------------------------------------------------------- */
     //! return physical vector of a cell of cubic pixels
     template <size_t dim>
     Eigen::Matrix<Real, dim, 1> get_vector(const Ccoord_t<dim> & ccoord, Real pix_size = 1.) {
       Eigen::Matrix<Real, dim, 1> retval;
       for (size_t i = 0; i < dim; ++i) {
         retval[i] = pix_size * ccoord[i];
       }
       return retval;
     }
 
     /* ---------------------------------------------------------------------- */
     //! return physical vector of a cell of general pixels
     template <size_t dim, typename T>
     Eigen::Matrix<T, dim, 1> get_vector(const Ccoord_t<dim> & ccoord,
                                         Eigen::Matrix<T, Dim_t(dim), 1> pix_size) {
       Eigen::Matrix<T, dim, 1> retval = pix_size;
       for (size_t i = 0; i < dim; ++i) {
         retval[i] *= ccoord[i];
       }
       return retval;
     }
 
     /* ---------------------------------------------------------------------- */
     //! return physical vector of a cell of general pixels
     template <size_t dim, typename T>
     Eigen::Matrix<T, dim, 1> get_vector(const Ccoord_t<dim> & ccoord,
                                         const std::array<T, dim>& pix_size) {
       Eigen::Matrix<T, dim, 1> retval{};
       for (size_t i = 0; i < dim; ++i) {
         retval[i] = pix_size[i]*ccoord[i];
       }
       return retval;
     }
 
 
     /* ---------------------------------------------------------------------- */
     //! get all strides from a row-major grid
     template <size_t dim>
     constexpr Ccoord_t<dim> get_default_strides(const Ccoord_t<dim> & sizes) {
       return internal::compute_strides<dim>(sizes,
                                             std::make_index_sequence<dim>{});
     }
 
     //----------------------------------------------------------------------------//
     //! get the i-th pixel in a grid of size sizes
     template <size_t dim>
     constexpr Ccoord_t<dim> get_ccoord(const Ccoord_t<dim> & resolutions,
                                        const Ccoord_t<dim> & locations,
-                                       Dim_t index, bool transposed=false) {
+                                       Dim_t index) {
       Ccoord_t<dim> retval{{0}};
       Dim_t factor{1};
       for (Dim_t i = dim-1; i >=0; --i) {
-        if (transposed) {
-          if (i == 0) {
-            retval[1] = index/factor%resolutions[0] + locations[0];
-          }
-          else if (i == 1) {
-            retval[0] = index/factor%resolutions[1] + locations[1];
-          }
-          else {
-            retval[i] = index/factor%resolutions[i] + locations[i];
-          }
-        }
-        else {
-          retval[i] = index/factor%resolutions[i] + locations[i];
-        }
+        retval[i] = index/factor%resolutions[i] + locations[i];
         if (i != 0 ) {
           factor *= resolutions[i];
         }
       }
       return retval;
     }
 
+    //----------------------------------------------------------------------------//
+    //! get the i-th pixel in a grid of size sizes
+    template <size_t dim, size_t... I>
+    constexpr Ccoord_t<dim> get_ccoord(const Ccoord_t<dim> & resolutions,
+                                       const Ccoord_t<dim> & locations,
+                                       Dim_t index,
+                                       std::index_sequence<I...>) {
+      Ccoord_t<dim> ccoord{get_ccoord<dim>(resolutions, locations, index)};
+      return Ccoord_t<dim>({ccoord[I]...});
+    }
+
     //----------------------------------------------------------------------------//
     //! get the linear index of a pixel in a given grid
     template <size_t dim>
     constexpr Dim_t get_index(const Ccoord_t<dim> & sizes,
                               const Ccoord_t<dim> & locations,
                               const Ccoord_t<dim> & ccoord) {
       Dim_t retval{0};
       Dim_t factor{1};
       for (Dim_t i = dim-1; i >=0; --i) {
         retval += (ccoord[i]-locations[i])*factor;
         if (i != 0) {
           factor *= sizes[i];
         }
       }
       return retval;
     }
 
     //----------------------------------------------------------------------------//
     //! get the linear index of a pixel given a set of strides
     template <size_t dim>
     constexpr Dim_t get_index_from_strides(const Ccoord_t<dim> & strides,
                                            const Ccoord_t<dim> & ccoord) {
       Dim_t retval{0};
       for (const auto & tup: akantu::zip(strides, ccoord)) {
         const auto & stride = std::get<0>(tup);
         const auto & ccord_ = std::get<1>(tup);
         retval += stride * ccord_;
       }
       return retval;
     }
 
     //----------------------------------------------------------------------------//
     //! get the number of pixels in a grid
     template <size_t dim>
     constexpr size_t get_size(const Ccoord_t<dim>& sizes) {
       Dim_t retval{1};
       for (size_t i = 0; i < dim; ++i) {
         retval *= sizes[i];
       }
       return retval;
     }
 
     //----------------------------------------------------------------------------//
     //! get the number of pixels in a grid given its strides
     template <size_t dim>
     constexpr size_t get_size_from_strides(const Ccoord_t<dim>& sizes,
                                            const Ccoord_t<dim>& strides) {
       return sizes[0]*strides[0];
     }
 
     /* ---------------------------------------------------------------------- */
     /**
      * centralises iterating over square (or cubic) discretisation grids
      */
-    template <size_t dim, bool transposed=false>
+    template <size_t dim, int ...dmap>
     class Pixels {
     public:
       //! constructor
       Pixels(const Ccoord_t<dim> & resolutions=Ccoord_t<dim>{},
              const Ccoord_t<dim> & locations=Ccoord_t<dim>{})
         :resolutions{resolutions}, locations{locations}{};
       //! copy constructor
       Pixels(const Pixels & other) = default;
       //! assignment operator
       Pixels & operator=(const Pixels & other) = default;
       virtual ~Pixels() = default;
 
       /**
        * iterators over `Pixels` dereferences to cell coordinates
        */
       class iterator
       {
       public:
         using value_type = Ccoord_t<dim>; //!< stl conformance
         using const_value_type = const value_type; //!< stl conformance
         using pointer = value_type*; //!< stl conformance
         using difference_type = std::ptrdiff_t; //!< stl conformance
         using iterator_category = std::forward_iterator_tag;//!<stl conformance
         using reference = value_type; //!< stl conformance
 
         //! constructor
         iterator(const Pixels & pixels, bool begin=true);
         virtual ~iterator() = default;
         //! dereferencing
         inline value_type operator*() const;
         //! pre-increment
         inline iterator & operator++();
         //! inequality
         inline bool operator!=(const iterator & other) const;
         //! equality
         inline bool operator==(const iterator & other) const;
 
       protected:
         const Pixels& pixels; //!< ref to pixels in cell
         size_t index; //!< index of currect pointed-to pixel
       };
       //! stl conformance
       inline iterator begin() const {return iterator(*this);}
       //! stl conformance
       inline iterator end() const {return iterator(*this, false);}
       //! stl conformance
       inline size_t size() const {return get_size(this->resolutions);}
     protected:
       Ccoord_t<dim> resolutions; //!< resolutions of this domain
       Ccoord_t<dim> locations; //!< locations of this domain
     };
 
     /* ---------------------------------------------------------------------- */
-    template <size_t dim, bool transposed>
-    Pixels<dim, transposed>::iterator::iterator(const Pixels & pixels, bool begin)
+    template <size_t dim, int ...dmap>
+    Pixels<dim, dmap...>::iterator::iterator(const Pixels & pixels, bool begin)
       :pixels{pixels}, index{begin? 0: get_size(pixels.resolutions)}
     {}
 
     /* ---------------------------------------------------------------------- */
-    template <size_t dim, bool transposed>
-    typename Pixels<dim, transposed>::iterator::value_type
-    Pixels<dim, transposed>::iterator::operator*() const {
+    template <size_t dim, int ...dmap>
+    typename Pixels<dim, dmap...>::iterator::value_type
+    Pixels<dim, dmap...>::iterator::operator*() const {
       return get_ccoord(pixels.resolutions, pixels.locations, this->index,
-                        transposed);
+                        std::conditional_t<sizeof...(dmap) == 0,
+                                           std::make_index_sequence<dim>,
+                                           std::index_sequence<dmap...>>{});
     }
 
     /* ---------------------------------------------------------------------- */
-    template <size_t dim, bool transposed>
+    template <size_t dim, int ...dmap>
     bool
-    Pixels<dim, transposed>::iterator::operator!=(const iterator &other) const {
+    Pixels<dim, dmap...>::iterator::operator!=(const iterator &other) const {
       return (this->index != other.index) || (&this->pixels != &other.pixels);
     }
 
     /* ---------------------------------------------------------------------- */
-    template <size_t dim, bool transposed>
+    template <size_t dim, int ...dmap>
     bool
-    Pixels<dim, transposed>::iterator::operator==(const iterator &other) const {
+    Pixels<dim, dmap...>::iterator::operator==(const iterator &other) const {
       return !(*this!= other);
     }
 
     /* ---------------------------------------------------------------------- */
-    template <size_t dim, bool transposed>
-    typename Pixels<dim, transposed>::iterator&
-    Pixels<dim, transposed>::iterator::operator++() {
+    template <size_t dim, int ...dmap>
+    typename Pixels<dim, dmap...>::iterator&
+    Pixels<dim, dmap...>::iterator::operator++() {
       ++this->index;
       return *this;
     }
 
   }  // CcoordOps
 
 }  // muSpectre
 
 #endif /* CCOORD_OPERATIONS_H */
diff --git a/src/fft/fftwmpi_engine.cc b/src/fft/fftwmpi_engine.cc
index 7d17f1f..1bf7c6f 100644
--- a/src/fft/fftwmpi_engine.cc
+++ b/src/fft/fftwmpi_engine.cc
@@ -1,223 +1,227 @@
 /**
  * @file   fftwmpi_engine.cc
  *
  * @author Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
  *
  * @date   06 Mar 2017
  *
  * @brief  implements the MPI-parallel fftw engine
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre 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
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "common/ccoord_operations.hh"
 #include "fft/fftwmpi_engine.hh"
 
 namespace muSpectre {
 
   template <Dim_t DimsS, Dim_t DimM>
   int FFTWMPIEngine<DimsS, DimM>::nb_engines{0};
 
   template <Dim_t DimS, Dim_t DimM>
   FFTWMPIEngine<DimS, DimM>::FFTWMPIEngine(Ccoord resolutions, Rcoord lengths,
                                            Communicator comm)
     :Parent{resolutions, lengths, comm}
   {
     if (!this->nb_engines) fftw_mpi_init();
     this->nb_engines++;
 
     std::array<ptrdiff_t, DimS> narr;
     std::copy(this->domain_resolutions.begin(), this->domain_resolutions.end(),
               narr.begin());
     narr[DimS-1] = this->domain_resolutions[DimS-1]/2+1;
     ptrdiff_t res_x, loc_x, res_y, loc_y;
     this->workspace_size =
       fftw_mpi_local_size_many_transposed(DimS, narr.data(),
                                           Field_t::nb_components,
                                           FFTW_MPI_DEFAULT_BLOCK,
                                           FFTW_MPI_DEFAULT_BLOCK,
                                           this->comm.get_mpi_comm(),
                                           &res_x, &loc_x, &res_y, &loc_y);
     this->fourier_resolutions[1] = this->fourier_resolutions[0];
     this->fourier_locations[1] = this->fourier_locations[0];
     this->resolutions[0] = res_x;
     this->locations[0] = loc_x;
     this->fourier_resolutions[0] = res_y;
     this->fourier_locations[0] = loc_y;
 
     for (auto & n: this->resolutions) {
       if (n == 0) {
         throw std::runtime_error("FFTW MPI planning returned zero resolution. "
                                  "You may need to run on fewer processes.");
       }
     }
     for (auto & n: this->fourier_resolutions) {
       if (n == 0) {
         throw std::runtime_error("FFTW MPI planning returned zero Fourier "
                                  "resolution. You may need to run on fewer "
                                  "processes.");
       }
     }
 
-    for (auto && pixel: CcoordOps::Pixels<DimS, true>(this->fourier_resolutions,
-                                                      this->fourier_locations)) {
+    for (auto && pixel:
+         std::conditional_t<
+           DimS==2,
+           CcoordOps::Pixels<DimS, 1, 0>,
+           CcoordOps::Pixels<DimS, 1, 0, 2>
+         >(this->fourier_resolutions, this->fourier_locations)) {
            this->work_space_container.add_pixel(pixel);
     }
   }
 
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void FFTWMPIEngine<DimS, DimM>::initialise(FFT_PlanFlags plan_flags) {
     if (this->initialised) {
       throw std::runtime_error("double initialisation, will leak memory");
     }
 
     // Initialize parent after local resolutions have been determined and
     // work space has been initialized
     Parent::initialise(plan_flags);
 
     this->real_workspace = fftw_alloc_real(2*this->workspace_size);
     // We need to check whether the workspace provided by our field is large
     // enough. MPI parallel FFTW may request a workspace size larger than the
     // nominal size of the complex buffer.
     if (long(this->work.size()*Field_t::nb_components) < this->workspace_size) {
       this->work.set_pad_size(this->workspace_size -
                               Field_t::nb_components*this->work.size());
     }
 
     unsigned int flags;
     switch (plan_flags) {
     case FFT_PlanFlags::estimate: {
       flags = FFTW_ESTIMATE;
       break;
     }
     case FFT_PlanFlags::measure: {
       flags = FFTW_MEASURE;
       break;
     }
     case FFT_PlanFlags::patient: {
       flags = FFTW_PATIENT;
       break;
     }
     default:
       throw std::runtime_error("unknown planner flag type");
       break;
     }
 
     std::array<ptrdiff_t, DimS> narr;
     std::copy(this->domain_resolutions.begin(), this->domain_resolutions.end(),
               narr.begin());
     Real * in = this->real_workspace;
     fftw_complex * out = reinterpret_cast<fftw_complex*>(this->work.data());
     this->plan_fft = fftw_mpi_plan_many_dft_r2c(
       DimS, narr.data(), Field_t::nb_components, FFTW_MPI_DEFAULT_BLOCK,
       FFTW_MPI_DEFAULT_BLOCK, in, out, this->comm.get_mpi_comm(),
       FFTW_MPI_TRANSPOSED_OUT | flags);
     if (this->plan_fft == nullptr) {
       throw std::runtime_error("r2c plan failed");
     }
 
     fftw_complex * i_in = reinterpret_cast<fftw_complex*>(this->work.data());
     Real * i_out = this->real_workspace;
 
     this->plan_ifft = fftw_mpi_plan_many_dft_c2r(
       DimS, narr.data(), Field_t::nb_components, FFTW_MPI_DEFAULT_BLOCK,
       FFTW_MPI_DEFAULT_BLOCK, i_in, i_out, this->comm.get_mpi_comm(),
       FFTW_MPI_TRANSPOSED_IN | flags);
     if (this->plan_ifft == nullptr) {
       throw std::runtime_error("c2r plan failed");
     }
     this->initialised = true;
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   FFTWMPIEngine<DimS, DimM>::~FFTWMPIEngine<DimS, DimM>() noexcept {
     if (this->real_workspace != nullptr) fftw_free(this->real_workspace);
     if (this->plan_fft != nullptr) fftw_destroy_plan(this->plan_fft);
     if (this->plan_ifft != nullptr) fftw_destroy_plan(this->plan_ifft);
     this->nb_engines--;
     if (!this->nb_engines) fftw_mpi_cleanup();
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   typename FFTWMPIEngine<DimS, DimM>::Workspace_t &
   FFTWMPIEngine<DimS, DimM>::fft (Field_t & field) {
     if (this->plan_fft == nullptr) {
       throw std::runtime_error("fft plan not initialised");
     }
     if (field.size() != CcoordOps::get_size(this->resolutions)) {
       throw std::runtime_error("size mismatch");
     }
     // Copy non-padded field to padded real_workspace.
     // Transposed output of M x N x L transform for >= 3 dimensions is padded
     // M x N x 2*(L/2+1).
     ptrdiff_t fstride = Field_t::nb_components*this->resolutions[DimS-1];
     ptrdiff_t wstride = Field_t::nb_components*2*(this->resolutions[DimS-1]/2+1);
     ptrdiff_t n = field.size()/this->resolutions[DimS-1];
 
     auto fdata = field.data();
     auto wdata = this->real_workspace;
     for (int i = 0; i < n; i++) {
       std::copy(fdata, fdata+fstride, wdata);
       fdata += fstride;
       wdata += wstride;
     }
     // Compute FFT
     fftw_mpi_execute_dft_r2c(
       this->plan_fft, this->real_workspace,
       reinterpret_cast<fftw_complex*>(this->work.data()));
     return this->work;
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void
   FFTWMPIEngine<DimS, DimM>::ifft (Field_t & field) const {
     if (this->plan_ifft == nullptr) {
       throw std::runtime_error("ifft plan not initialised");
     }
     if (field.size() != CcoordOps::get_size(this->resolutions)) {
       throw std::runtime_error("size mismatch");
     }
     // Compute inverse FFT
     fftw_mpi_execute_dft_c2r(
       this->plan_ifft, reinterpret_cast<fftw_complex*>(this->work.data()),
       this->real_workspace);
     // Copy non-padded field to padded real_workspace.
     // Transposed output of M x N x L transform for >= 3 dimensions is padded
     // M x N x 2*(L/2+1).
     ptrdiff_t fstride = Field_t::nb_components*this->resolutions[DimS-1];
     ptrdiff_t wstride = Field_t::nb_components*2*(this->resolutions[DimS-1]/2+1);
     ptrdiff_t n = field.size()/this->resolutions[DimS-1];
 
     auto fdata = field.data();
     auto wdata = this->real_workspace;
     for (int i = 0; i < n; i++) {
       std::copy(wdata, wdata+fstride, fdata);
       fdata += fstride;
       wdata += wstride;
     }
   }
 
   template class FFTWMPIEngine<twoD, twoD>;
   template class FFTWMPIEngine<twoD, threeD>;
   template class FFTWMPIEngine<threeD, threeD>;
 }  // muSpectre
diff --git a/src/fft/pfft_engine.cc b/src/fft/pfft_engine.cc
index a810d08..b827223 100644
--- a/src/fft/pfft_engine.cc
+++ b/src/fft/pfft_engine.cc
@@ -1,235 +1,242 @@
 /**
  * @file   pfft_engine.cc
  *
  * @author Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
  *
  * @date   06 Mar 2017
  *
  * @brief  implements the MPI-parallel pfft engine
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre 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
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "common/ccoord_operations.hh"
 #include "fft/pfft_engine.hh"
 
 namespace muSpectre {
 
   template <Dim_t DimsS, Dim_t DimM>
   int PFFTEngine<DimsS, DimM>::nb_engines{0};
 
   template <Dim_t DimS, Dim_t DimM>
   PFFTEngine<DimS, DimM>::PFFTEngine(Ccoord resolutions, Rcoord lengths,
                                      Communicator comm)
     :Parent{resolutions, lengths, comm}
   {
     if (!this->nb_engines) pfft_init();
     this->nb_engines++;
 
     int size{comm.size()};
     int dimx{size};
     int dimy{1};
-    if (DimS > 2) {
+    //if (DimS > 2) {
+    if (false) {
       dimy = int(sqrt(size));
       while ((size/dimy)*dimy != size)  dimy--;
       dimx = size/dimy;
     }
 
     std::cout << "PFFT process mesh: " << dimx << "x" << dimy << std::endl;
 
-    if (DimS > 2) {
+    //if (DimS > 2) {
+    if (false) {
       if (pfft_create_procmesh_2d(this->comm.get_mpi_comm(), dimx, dimy,
                                   &this->mpi_comm)) {
         throw std::runtime_error("Failed to create 2d PFFT process mesh.");
       }
     }
     else {
       this->mpi_comm = this->comm.get_mpi_comm();
     }
 
     // FIXME! Code presently always use a one-dimensional process mesh but
     // should use a two-dimensional process mesh for three-dimensional data.
     std::array<ptrdiff_t, DimS> narr;
     std::copy(this->domain_resolutions.begin(), this->domain_resolutions.end(),
               narr.begin());
     ptrdiff_t res[DimS], loc[DimS], fres[DimS], floc[DimS];
     this->workspace_size =
       pfft_local_size_many_dft_r2c(DimS, narr.data(), narr.data(), narr.data(),
                                    Field_t::nb_components,
                                    PFFT_DEFAULT_BLOCK, PFFT_DEFAULT_BLOCKS,
                                    this->mpi_comm,
                                    PFFT_TRANSPOSED_OUT,
                                    res, loc, fres, floc);
     std::copy(res, res+DimS, this->resolutions.begin());
     std::copy(loc, loc+DimS, this->locations.begin());
     std::copy(fres, fres+DimS, this->fourier_resolutions.begin());
     std::copy(floc, floc+DimS, this->fourier_locations.begin());
-    for (int i = 0; i < DimS-1; i++) {
+    //for (int i = 0; i < DimS-1; i++) {
+    for (int i = 0; i < 1; i++) {
       std::swap(this->fourier_resolutions[i], this->fourier_resolutions[i+1]);
       std::swap(this->fourier_locations[i], this->fourier_locations[i+1]);
     }
 
     for (auto & n: this->resolutions) {
       if (n == 0) {
         throw std::runtime_error("PFFT planning returned zero resolution. "
                                  "You may need to run on fewer processes.");
       }
     }
     for (auto & n: this->fourier_resolutions) {
       if (n == 0) {
         throw std::runtime_error("PFFT planning returned zero Fourier "
                                  "resolution. You may need to run on fewer "
                                  "processes.");
       }
     }
 
-    for (auto && pixel: CcoordOps::Pixels<DimS, true>(this->fourier_resolutions,
-                                                      this->fourier_locations)) {
+    for (auto && pixel:
+         std::conditional_t<
+         DimS==2,
+         CcoordOps::Pixels<DimS, 1, 0>,
+         CcoordOps::Pixels<DimS, 1, 0, 2>
+         >(this->fourier_resolutions, this->fourier_locations)) {
       this->work_space_container.add_pixel(pixel);
     }
   }
 
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void PFFTEngine<DimS, DimM>::initialise(FFT_PlanFlags plan_flags) {
     if (this->initialised) {
       throw std::runtime_error("double initialisation, will leak memory");
     }
 
     // Initialize parent after local resolutions have been determined and
     // work space has been initialized
     Parent::initialise(plan_flags);
 
     this->real_workspace = pfft_alloc_real(2*this->workspace_size);
     // We need to check whether the workspace provided by our field is large
     // enough. PFFT may request a workspace size larger than the nominal size
     // of the complex buffer.
     if (long(this->work.size()*Field_t::nb_components) < this->workspace_size) {
       this->work.set_pad_size(this->workspace_size -
                               Field_t::nb_components*this->work.size());
     }
 
     unsigned int flags;
     switch (plan_flags) {
     case FFT_PlanFlags::estimate: {
       flags = PFFT_ESTIMATE;
       break;
     }
     case FFT_PlanFlags::measure: {
       flags = PFFT_MEASURE;
       break;
     }
     case FFT_PlanFlags::patient: {
       flags = PFFT_PATIENT;
       break;
     }
     default:
       throw std::runtime_error("unknown planner flag type");
       break;
     }
 
     std::array<ptrdiff_t, DimS> narr;
     std::copy(this->domain_resolutions.begin(), this->domain_resolutions.end(),
               narr.begin());
     Real * in = this->real_workspace;
     pfft_complex * out = reinterpret_cast<pfft_complex*>(this->work.data());
     this->plan_fft = pfft_plan_many_dft_r2c(DimS, narr.data(), narr.data(),
                                             narr.data(), Field_t::nb_components,
                                             PFFT_DEFAULT_BLOCKS,
                                             PFFT_DEFAULT_BLOCKS,
                                             in, out, this->mpi_comm,
                                             PFFT_FORWARD,
                                             PFFT_TRANSPOSED_OUT | flags);
     if (this->plan_fft == nullptr) {
       throw std::runtime_error("r2c plan failed");
     }
 
     pfft_complex * i_in = reinterpret_cast<pfft_complex*>(this->work.data());
     Real * i_out = this->real_workspace;
 
     this->plan_ifft = pfft_plan_many_dft_c2r(DimS, narr.data(), narr.data(),
                                              narr.data(),
                                              Field_t::nb_components,
                                              PFFT_DEFAULT_BLOCKS,
                                              PFFT_DEFAULT_BLOCKS,
                                              i_in, i_out,
                                              this->mpi_comm,
                                              PFFT_BACKWARD,
                                              PFFT_TRANSPOSED_IN | flags);
     if (this->plan_ifft == nullptr) {
       throw std::runtime_error("c2r plan failed");
     }
     this->initialised = true;
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   PFFTEngine<DimS, DimM>::~PFFTEngine<DimS, DimM>() noexcept {
     if (this->real_workspace != nullptr) pfft_free(this->real_workspace);
     if (this->plan_fft != nullptr) pfft_destroy_plan(this->plan_fft);
     if (this->plan_ifft != nullptr) pfft_destroy_plan(this->plan_ifft);
     if (this->mpi_comm != this->comm.get_mpi_comm()) {
       MPI_Comm_free(&this->mpi_comm);
     }
     this->nb_engines--;
     if (!this->nb_engines) pfft_cleanup();
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   typename PFFTEngine<DimS, DimM>::Workspace_t &
   PFFTEngine<DimS, DimM>::fft (Field_t & field) {
     if (!this->plan_fft) {
         throw std::runtime_error("fft plan not allocated");
     }
     if (field.size() != CcoordOps::get_size(this->resolutions)) {
       throw std::runtime_error("size mismatch");
     }
     // Copy field data to workspace buffer. This is necessary because workspace
     // buffer is larger than field buffer.
     std::copy(field.data(), field.data()+Field_t::nb_components*field.size(),
               this->real_workspace);
     pfft_execute_dft_r2c(
       this->plan_fft, this->real_workspace,
       reinterpret_cast<pfft_complex*>(this->work.data()));
     return this->work;
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void
   PFFTEngine<DimS, DimM>::ifft (Field_t & field) const {
     if (!this->plan_ifft) {
         throw std::runtime_error("ifft plan not allocated");
     }
     if (field.size() != CcoordOps::get_size(this->resolutions)) {
       throw std::runtime_error("size mismatch");
     }
     pfft_execute_dft_c2r(
       this->plan_ifft, reinterpret_cast<pfft_complex*>(this->work.data()),
       this->real_workspace);
     std::copy(this->real_workspace,
               this->real_workspace+Field_t::nb_components*field.size(),
               field.data());
   }
 
   template class PFFTEngine<twoD, twoD>;
   template class PFFTEngine<twoD, threeD>;
   template class PFFTEngine<threeD, threeD>;
 }  // muSpectre