diff --git a/.gitignore b/.gitignore
index a616d26..e15f15a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
 /build/CMake*
 *~
 /build/*
+/build-o3/
diff --git a/src/common/ccoord_operations.hh b/src/common/ccoord_operations.hh
index f46911b..16d9d09 100644
--- a/src/common/ccoord_operations.hh
+++ b/src/common/ccoord_operations.hh
@@ -1,84 +1,84 @@
 /**
  * file   ccoord_operations.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   29 Sep 2017
  *
  * @brief  common operations on pixel addressing 
  *
  * @section LICENCE
  *
  * Copyright (C) 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 "common/common.hh"
 
 #ifndef CCOORD_OPERATIONS_H
 #define CCOORD_OPERATIONS_H
 
 namespace muSpectre {
 
   namespace CcoordOps {
 
     //----------------------------------------------------------------------------//
     template <size_t dim>
     constexpr Ccoord_t<dim> get_ccoord(const Ccoord_t<dim> & sizes, Dim_t index) {
       Ccoord_t<dim> retval{0};
       Dim_t factor{1};
       for (Dim_t i = dim-1; i >=0; --i) {
         retval[i] = index/factor%sizes[i];
-        if (dim != 0 ) {
+        if (i != 0 ) {
           factor *= sizes[i];
         }
       }
       return retval;
     }
 
     //----------------------------------------------------------------------------//
     template <size_t dim>
     constexpr Dim_t get_index(const Ccoord_t<dim> & sizes,
                               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]*factor;
-        if (dim != 0) {
+        if (i != 0) {
           factor *= sizes[i];
         }
       }
       return retval;
     }
 
     //----------------------------------------------------------------------------//
     template <size_t dim>
     constexpr size_t get_size(const Ccoord_t<dim>& sizes) {
       return std::accumulate(sizes.begin(), sizes.end(), 1,
                              std::multiplies<size_t>());
     }
 
 
 
   }  // CcoordOps
 
 }  // muSpectre
 
 #endif /* CCOORD_OPERATIONS_H */
diff --git a/src/common/field.hh b/src/common/field.hh
index 97b9b11..c443b8e 100644
--- a/src/common/field.hh
+++ b/src/common/field.hh
@@ -1,408 +1,408 @@
 /**
  * file   field.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   07 Sep 2017
  *
  * @brief  header-only implementation of a field for field collections
  *
  * @section LICENCE
  *
  * Copyright (C) 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.
  */
 
 
 #ifndef FIELD_H
 #define FIELD_H
 
 #include <string>
 #include <utility>
 #include <typeinfo>
 #include <vector>
 #include <algorithm>
 #include <Eigen/Dense>
 #include <cmath>
 #include <memory>
 #include <type_traits>
 
 namespace muSpectre {
 
   namespace internal {
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection>
     class FieldBase
     {
 
     protected:
       //! constructor
       //! unique name (whithin Collection)
       //! number of components
       //! collection to which this field belongs (eg, material, system)
       FieldBase(std::string unique_name,
                 size_t nb_components,
                 FieldCollection & collection);
 
     public:
 
       //! Copy constructor
       FieldBase(const FieldBase &other) = delete;
 
       //! Move constructor
       FieldBase(FieldBase &&other) noexcept = delete;
 
       //! Destructor
       virtual ~FieldBase() noexcept = default;
 
       //! Copy assignment operator
       FieldBase& operator=(const FieldBase &other) = delete;
 
       //! Move assignment operator
       FieldBase& operator=(FieldBase &&other) noexcept = delete;
 
       /* ---------------------------------------------------------------------- */
       //!Identifying accessors
       //! return field name
       inline const std::string & get_name() const;
       //! return field type
       //inline const Field_t & get_type() const;
       //! return my collection (for iterating)
       inline const FieldCollection & get_collection() const;
       //! return my collection (for iterating)
       inline const size_t & get_nb_components() const;
       //! return type_id of stored type
       virtual const std::type_info & get_stored_typeid() const = 0;
 
       virtual size_t size() const = 0;
 
       //! initialise field to zero (do more complicated initialisations through
       //! fully typed maps)
       virtual void set_zero() = 0;
 
       //! give access to collections
       friend FieldCollection;
 
     protected:
       /* ---------------------------------------------------------------------- */
       //! allocate memory etc
       virtual void resize(size_t size) = 0;
       const std::string name;
       const size_t nb_components;
       const FieldCollection & collection;
     private:
     };
 
 
     /* ---------------------------------------------------------------------- */
     //! declaraton for friending
     template <class FieldCollection, typename T, Dim_t NbComponents>
     class FieldMap;
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     class TypedFieldBase: public FieldBase<FieldCollection>
     {
       friend class FieldMap<FieldCollection, T, NbComponents>;
     public:
       using Parent = FieldBase<FieldCollection>;
       using Base = Parent;
       //using storage_type = Eigen::Array<T, Eigen::Dynamic, NbComponents>;
       using StoredType = Eigen::Array<T, NbComponents, 1>;
       using StorageType = std::vector<StoredType,
                                        Eigen::aligned_allocator<StoredType>>;
       TypedFieldBase(std::string unique_name,
                      FieldCollection& collection);
       virtual ~TypedFieldBase() = default;
       //! return type_id of stored type
       virtual const std::type_info & get_stored_typeid() const;
 
       //! initialise field to zero (do more complicated initialisations through
       //! fully typed maps)
       inline void set_zero() override final;
       inline void push_back(const StoredType & value);
       template< class... Args>
       inline void emplace_back(Args&&... args);
       size_t size() const override final;
     protected:
       inline T* get_ptr_to_entry(const size_t&& index);
       inline T& get_ref_to_entry(const size_t&& index);
       inline virtual void resize(size_t size) override final;
       StorageType array;
     };
 
   }  // internal
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   class TensorField: public internal::TypedFieldBase<FieldCollection,
                                                      T,
                                                      ipow(dim,order)>
   {
   public:
     using Parent = internal::TypedFieldBase<FieldCollection,
                                             T,
                                             ipow(dim,order)>;
     using Base = typename Parent::Base;
-    using Field_p = std::unique_ptr<internal::FieldBase<FieldCollection>>;
+    using Field_p = typename FieldCollection::Field_p;
     using component_type = T;
     //! Copy constructor
     TensorField(const TensorField &other) = delete;
 
     //! Move constructor
     TensorField(TensorField &&other) noexcept = delete;
 
     //! Destructor
     virtual ~TensorField() noexcept = default;
 
     //! Copy assignment operator
     TensorField& operator=(const TensorField &other) = delete;
 
     //! Move assignment operator
     TensorField& operator=(TensorField &&other) noexcept = delete;
 
     //! accessors
     inline Dim_t get_order() const;
     inline Dim_t get_dim() const;
 
 
     //! factory function
     template<class FieldType, class CollectionType, typename... Args>
     friend typename FieldType::Base&  make_field(std::string unique_name,
                              CollectionType & collection,
                              Args&&... args);
 
   protected:
     //! constructor protected!
     TensorField(std::string unique_name, FieldCollection & collection);
 
   private:
   };
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol=NbRow>
   class MatrixField: public internal::TypedFieldBase<FieldCollection,
                                                      T,
                                                      NbRow*NbCol>
   {
   public:
     using Parent = internal::TypedFieldBase<FieldCollection,
                                             T,
                                             NbRow*NbCol>;
     using Base = typename Parent::Base;
     using Field_p = std::unique_ptr<internal::FieldBase<FieldCollection>>;
     using component_type = T;
     //! Copy constructor
     MatrixField(const MatrixField &other) = delete;
 
     //! Move constructor
     MatrixField(MatrixField &&other) noexcept = delete;
 
     //! Destructor
     virtual ~MatrixField() noexcept = default;
 
     //! Copy assignment operator
     MatrixField& operator=(const MatrixField &other) = delete;
 
     //! Move assignment operator
     MatrixField& operator=(MatrixField &&other) noexcept = delete;
 
     //! accessors
     inline Dim_t get_nb_row() const;
     inline Dim_t get_nb_col() const;
 
 
     //! factory function
     template<class FieldType, class CollectionType, typename... Args>
     friend typename FieldType::Base&  make_field(std::string unique_name,
                              CollectionType & collection,
                              Args&&... args);
 
   protected:
     //! constructor protected!
     MatrixField(std::string unique_name, FieldCollection & collection);
 
   private:
   };
 
   /* ---------------------------------------------------------------------- */
   //! convenience alias
   template <class FieldCollection, typename T>
   using ScalarField = TensorField<FieldCollection, T, 1, 1>;
   /* ---------------------------------------------------------------------- */
   // Implementations
   namespace internal {
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection>
     FieldBase<FieldCollection>::FieldBase(std::string unique_name,
                                           size_t nb_components_,
                                           FieldCollection & collection_)
       :name(unique_name), nb_components(nb_components_),
        collection(collection_) {}
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection>
     inline const std::string & FieldBase<FieldCollection>::get_name() const {
       return this->name;
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection>
     inline const FieldCollection & FieldBase<FieldCollection>::
     get_collection() const {
       return this->collection;
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection>
     inline const size_t & FieldBase<FieldCollection>::
     get_nb_components() const {
       return this->nb_components;
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     TypedFieldBase<FieldCollection, T, NbComponents>::
     TypedFieldBase(std::string unique_name, FieldCollection & collection)
       :FieldBase<FieldCollection>(unique_name, NbComponents, collection){
       static_assert
         ((std::is_arithmetic<T>::value ||
           std::is_same<T, Complex>::value),
          "Use TypedFieldBase for integer, real or complex scalars for T");
       static_assert(NbComponents > 0,
                     "Only fields with more than 0 components");
     }
 
     /* ---------------------------------------------------------------------- */
     //! return type_id of stored type
     template <class FieldCollection, typename T, Dim_t NbComponents>
     const std::type_info & TypedFieldBase<FieldCollection, T, NbComponents>::
     get_stored_typeid() const {
       return typeid(T);
     }
 
     /* ---------------------------------------------------------------------- */
     template<class FieldCollection, typename T, Dim_t NbComponents>
     void
     TypedFieldBase<FieldCollection, T, NbComponents>::
     set_zero() {
       std::fill(this->array.begin(), this->array.end(), T{});
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     size_t TypedFieldBase<FieldCollection, T, NbComponents>::
     size() const {
       return this->array.size();
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     T* TypedFieldBase<FieldCollection, T, NbComponents>::
     get_ptr_to_entry(const size_t&& index) {
       return &this->array[std::move(index)](0, 0);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     T& TypedFieldBase<FieldCollection, T, NbComponents>::
     get_ref_to_entry(const size_t && index) {
       return this->array[std::move(index)](0, 0);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     void TypedFieldBase<FieldCollection, T, NbComponents>::
     resize(size_t size) {
       this->array.resize(size);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     void TypedFieldBase<FieldCollection, T, NbComponents>::
     push_back(const StoredType & value) {
       this->array.push_back(value);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     template <class... Args>
     void TypedFieldBase<FieldCollection, T, NbComponents>::
     emplace_back(Args&&... args) {
       this->array.emplace_back(std::move(args...));
     }
 
   }  // internal
 
   /* ---------------------------------------------------------------------- */
   //! Factory function, guarantees that only fields get created that are
   //! properly registered and linked to a collection.
   template<class FieldType, class FieldCollection, typename... Args>
   typename FieldType::Base &
   make_field(std::string unique_name,
              FieldCollection & collection,
              Args&&... args) {
     auto && ptr = std::unique_ptr<FieldType>{
       new FieldType(unique_name, collection, args...)};
     auto& retref = *ptr;
     collection.register_field(std::move(ptr));
     return retref;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   TensorField<FieldCollection, T, order, dim>::
   TensorField(std::string unique_name, FieldCollection & collection)
     :Parent(unique_name, collection) {}
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   Dim_t TensorField<FieldCollection, T, order, dim>::
   get_order() const {
     return order;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   Dim_t TensorField<FieldCollection, T, order, dim>::
   get_dim() const {
     return dim;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   MatrixField<FieldCollection, T, NbRow, NbCol>::
   MatrixField(std::string unique_name, FieldCollection & collection)
     :Parent(unique_name, collection) {}
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   Dim_t MatrixField<FieldCollection, T, NbRow, NbCol>::
   get_nb_col() const {
     return NbCol;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   Dim_t MatrixField<FieldCollection, T, NbRow, NbCol>::
   get_nb_row() const {
     return NbRow;
   }
 
 }  // muSpectre
 #endif /* FIELD_H */
diff --git a/src/common/field_collection.hh b/src/common/field_collection.hh
index c002656..8842f73 100644
--- a/src/common/field_collection.hh
+++ b/src/common/field_collection.hh
@@ -1,307 +1,48 @@
 /**
  * file   field_collection.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   07 Sep 2017
  *
  * @brief  Provides pixel-iterable containers for scalar and tensorial fields,
  *         addressable by field name
  *
  * @section LICENCE
  *
  * Copyright (C) 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.
  */
 
 #ifndef FIELD_COLLECTION_H
 #define FIELD_COLLECTION_H
 
-#include "common/common.hh"
-#include "common/ccoord_operations.hh"
-#include "common/field.hh"
-#include <memory>
-#include <exception>
-#include <string>
-#include <sstream>
-#include <map>
-#include <vector>
-#include <algorithm>
-#include <numeric>
-#include <iostream>
+#include "common/field_collection_global.hh"
+#include "common/field_collection_local.hh"
 
 namespace muSpectre {
 
-  class FieldCollectionError: public std::runtime_error {
-  public:
-    explicit FieldCollectionError(const std::string& what)
-      :std::runtime_error(what){}
-    explicit FieldCollectionError(const char * what)
-      :std::runtime_error(what){}
-  };
-
-  class FieldError: public FieldCollectionError {
-    using Parent = FieldCollectionError;
-  public:
-    explicit FieldError(const std::string& what)
-      :Parent(what){}
-    explicit FieldError(const char * what)
-      :Parent(what){}
-  };
-  class FieldInterpretationError: public FieldError
-  {
-  public:
-    explicit FieldInterpretationError(const std::string & what)
-      :FieldError(what){}
-    explicit FieldInterpretationError(const char * what)
-      :FieldError(what){}
-  };
-
-
-
-  /* ---------------------------------------------------------------------- */
-  //! DimS spatial dimension (dimension of problem
-  //! DimM material_dimension (dimension of constitutive law)
-  //! Global determines whether this field is present everywhere or per material
   template <Dim_t DimS, Dim_t DimM, bool Global = true>
-  class FieldCollection
-  {
-  public:
-    using Field = internal::FieldBase<FieldCollection>;
-    using Field_p = std::unique_ptr<Field>;
-    using Ccoord = Ccoord_t<DimS>;
-
-    //! Default constructor
-    FieldCollection();
-
-    //! Copy constructor
-    FieldCollection(const FieldCollection &other) = delete;
-
-    //! Move constructor
-    FieldCollection(FieldCollection &&other) noexcept = delete;
-
-    //! Destructor
-    virtual ~FieldCollection() noexcept = default;
-
-    //! Copy assignment operator
-    FieldCollection& operator=(const FieldCollection &other) = delete;
-
-    //! Move assignment operator
-    FieldCollection& operator=(FieldCollection &&other) noexcept = delete;
-
-    //! add a pixel/voxel to the field collection
-    template <bool NotGlobal = !Global>
-    inline std::enable_if_t<NotGlobal> add_pixel(const Ccoord & local_ccoord);
-
-    //! allocate memory, etc. At this point, the field knows how many
-    //! entries it should have (global fields know from the 'sizes'
-    //! parameter and non-global fields know from the size of the
-    //! pixel container, which was grown using 'add_pixel(...)'. The
-    //! job of initialise is to check that all fields are either of
-    //! size 0, in which case they need to be allocated, or are of the
-    //! same size as the pixel_container. Any field of a different
-    //! size is wrong. TODO: check whether it makes sense to put a
-    //! runtime check here
-    template <bool NotGlobal = !Global>
-    inline std::enable_if_t<NotGlobal> initialise();
-    template <bool isGlobal = Global>
-    inline std::enable_if_t<isGlobal> initialise(Ccoord sizes);
-
-    //! Register a new field (fields need to be in heap, so I want to keep them
-    //! yas shared pointers
-    void register_field(Field_p&& field);
-
-    //! for return values of iterators
-    constexpr inline static Dim_t spatial_dim();
-
-    //! for return values of iterators
-    inline Dim_t get_spatial_dim() const;
-
-    //! retrieve field by unique_name
-    inline Field& operator[](std::string unique_name);
-
-    //! retrieve field by unique_name with bounds checking
-    inline Field& at(std::string unique_name);
-
-    //! returns size of collection, this refers to the number of pixels handled
-    //! by the collection, not the number of fields
-    inline size_t size() const {return this->size_;}
-
-    //! return the pixel sizes
-    template <bool isGlobal = Global>
-    inline const std::enable_if_t<isGlobal, Ccoord> & get_sizes() const;
-
-    //! returns the linear index corresponding to cell coordinates
-    inline size_t get_index(Ccoord && ccoord) const;
-    //! returns the cell coordinates corresponding to a linear index
-    inline Ccoord get_ccoord(size_t index) const;
-
-  protected:
-    std::map<const std::string, Field_p> fields;
-    bool is_initialised = false;
-    const Uint id;
-    static Uint counter;
-    size_t size_ = 0;
-    Ccoord sizes;
-    //! container of pixel coords for non-global collections
-    std::vector<Ccoord> ccoords;
-    //! container of indices for non-global collections (slow!)
-    std::map<Ccoord, size_t> indices;
-  private:
-  };
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  Uint FieldCollection<DimS, DimM, Global>::counter{0};
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  FieldCollection<DimS, DimM, Global>::FieldCollection()
-    :id(counter++){}
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  template <bool NotGlobal>
-  std::enable_if_t<NotGlobal> FieldCollection<DimS, DimM, Global>::
-  add_pixel(const Ccoord & local_ccoord) {
-    this->indices[local_ccoord] = this->ccoords.size();
-    this->ccoords.push_back(local_ccoord);
-  }
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  template <bool NotGlobal>
-  std::enable_if_t<NotGlobal> FieldCollection<DimS, DimM, Global>::
-  initialise() {
-    this->size_ = this->ccoords.size();
-    std::for_each(std::begin(this->fields), std::end(this->fields),
-                  [this](auto && item) {
-                    auto && field = *item.second;
-                    const auto field_size = field.size();
-                    if (field_size == 0) {
-                      field.resize(this->size());
-                    } else if (field_size != this->size()) {
-                      std::stringstream err_stream;
-                      err_stream << "Field '" << field.get_name()
-                                 << "' contains " << field_size
-                                 << " entries, but the field collection "
-                                 << "has " << this->size() << " pixels";
-                      throw FieldCollectionError(err_stream.str());
-                    }
-                  });
-    this->is_initialised = true;
-  }
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  template <bool isGlobal>
-  std::enable_if_t<isGlobal> FieldCollection<DimS, DimM, Global>::
-  initialise(Ccoord sizes) {
-    this->size_ = std::accumulate(sizes.begin(), sizes.end(), 1,
-                                   std::multiplies<Dim_t>());
-    this->sizes = sizes;
-    std::for_each(std::begin(this->fields), std::end(this->fields),
-                  [this](auto && item) {
-                    auto && field = *item.second;
-                    const auto field_size = field.size();
-                    if (field_size == 0) {
-                      field.resize(this->size());
-                    } else if (field_size != this->size()) {
-                      std::stringstream err_stream;
-                      err_stream << "Field '" << field.get_name()
-                                 << "' contains " << field_size
-                                 << " entries, but the field collection "
-                                 << "has " << this->size() << " pixels";
-                      throw FieldCollectionError(err_stream.str());
-                    }
-                  });
-    this->is_initialised = true;
-  }
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  void FieldCollection<DimS, DimM, Global>::register_field(Field_p &&field) {
-    auto&& search_it = this->fields.find(field->get_name());
-    auto&& does_exist = search_it != this->fields.end();
-    if (does_exist) {
-      std::stringstream err_str;
-      err_str << "a field named " << field->get_name()
-              << "is already registered in this field collection. "
-              << "Currently registered fields: ";
-      for (const auto& name_field_pair: this->fields) {
-        err_str << ", " << name_field_pair.first;
-      }
-      throw FieldCollectionError(err_str.str());
-    }
-    this->fields[field->get_name()] = std::move(field);
-  }
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  constexpr Dim_t FieldCollection<DimS, DimM, Global>::spatial_dim() {
-    return DimS;
-  }
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  Dim_t FieldCollection<DimS, DimM, Global>::get_spatial_dim() const {
-    return DimS;
-  }
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  typename FieldCollection<DimS, DimM, Global>::Field&
-  FieldCollection<DimS, DimM, Global>::operator[](std::string unique_name) {
-    return *(this->fields[unique_name]);
-  }
-
-  /* ---------------------------------------------------------------------- */
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  typename FieldCollection<DimS, DimM, Global>::Field&
-  FieldCollection<DimS, DimM, Global>::at(std::string unique_name) {
-    return *(this->fields.at(unique_name));
-  }
-
-  //----------------------------------------------------------------------------//
-  //! return the pixel sizes
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  template <bool isGlobal>
-  const std::enable_if_t<isGlobal, Ccoord_t<DimS>>&
-  FieldCollection<DimS, DimM, Global>::get_sizes() const {
-    return this->sizes;
-  }
-
-  //----------------------------------------------------------------------------//
-  //! returns the linear index corresponding to cell coordinates
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  size_t
-  FieldCollection<DimS, DimM, Global>::get_index(Ccoord && ccoord) const {
-    return CcoordOps::get_index(this->sizes, std::move(ccoord));
-  }
-
-  //----------------------------------------------------------------------------//
-  //! returns the cell coordinates corresponding to a linear index
-  template <Dim_t DimS, Dim_t DimM, bool Global>
-  typename FieldCollection<DimS, DimM, Global>::Ccoord
-  FieldCollection<DimS, DimM, Global>::get_ccoord(size_t index) const {
-    return CcoordOps::get_ccoord(this->sizes, std::move(index));
-  }
+  using FieldCollection = std::conditional_t<
+    Global,
+    GlobalFieldCollection<DimS, DimM>,
+    LocalFieldCollection<DimS, DimM>>;
 
 
 }  // muSpectre
 
 #endif /* FIELD_COLLECTION_H */
diff --git a/src/common/field_collection_base.hh b/src/common/field_collection_base.hh
new file mode 100644
index 0000000..9202036
--- /dev/null
+++ b/src/common/field_collection_base.hh
@@ -0,0 +1,185 @@
+/**
+ * file   field_collection_base.hh
+ *
+ * @author Till Junge <till.junge@altermail.ch>
+ *
+ * @date   05 Nov 2017
+ *
+ * @brief  Base class for field collections
+ *
+ * @section LICENCE
+ *
+ * Copyright (C) 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.
+ */
+
+#ifndef FIELD_COLLECTION_BASE_H
+#define FIELD_COLLECTION_BASE_H
+
+#include <map>
+#include <vector>
+#include "common/common.hh"
+#include "common/field.hh"
+
+
+namespace muSpectre {
+  /* ---------------------------------------------------------------------- */
+  class FieldCollectionError: public std::runtime_error {
+  public:
+    explicit FieldCollectionError(const std::string& what)
+      :std::runtime_error(what){}
+    explicit FieldCollectionError(const char * what)
+      :std::runtime_error(what){}
+  };
+
+  class FieldError: public FieldCollectionError {
+    using Parent = FieldCollectionError;
+  public:
+    explicit FieldError(const std::string& what)
+      :Parent(what){}
+    explicit FieldError(const char * what)
+      :Parent(what){}
+  };
+  class FieldInterpretationError: public FieldError
+  {
+  public:
+    explicit FieldInterpretationError(const std::string & what)
+      :FieldError(what){}
+    explicit FieldInterpretationError(const char * what)
+      :FieldError(what){}
+  };
+
+  /* ---------------------------------------------------------------------- */
+  //! DimS spatial dimension (dimension of problem)
+  //! DimM material_dimension (dimension of constitutive law)
+  template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
+  class FieldCollectionBase
+  {
+  public:
+    using Field = internal::FieldBase<FieldCollectionDerived>;
+    using Field_p = std::unique_ptr<Field>;
+    using Ccoord = Ccoord_t<DimS>;
+
+    //! Default constructor
+    FieldCollectionBase();
+
+    //! Copy constructor
+    FieldCollectionBase(const FieldCollectionBase &other) = delete;
+
+    //! Move constructor
+    FieldCollectionBase(FieldCollectionBase &&other) noexcept = delete;
+
+    //! Destructor
+    virtual ~FieldCollectionBase() noexcept = default;
+
+    //! Copy assignment operator
+    FieldCollectionBase& operator=(const FieldCollectionBase &other) = delete;
+
+    //! Move assignment operator
+    FieldCollectionBase& operator=(FieldCollectionBase &&other) noexcept = delete;
+
+    //! Register a new field (fields need to be in heap, so I want to keep them
+    //! yas shared pointers
+    void register_field(Field_p&& field);
+
+    //! for return values of iterators
+    constexpr inline static Dim_t spatial_dim();
+
+    //! for return values of iterators
+    inline Dim_t get_spatial_dim() const;
+
+    //! retrieve field by unique_name
+    inline Field& operator[](std::string unique_name);
+
+    //! retrieve field by unique_name with bounds checking
+    inline Field& at(std::string unique_name);
+
+    //! returns size of collection, this refers to the number of pixels handled
+    //! by the collection, not the number of fields
+    inline size_t size() const {return this->size_;}
+
+  protected:
+    std::map<const std::string, Field_p> fields;
+    bool is_initialised = false;
+    const Uint id;
+    static Uint counter;
+    size_t size_ = 0;
+  private:
+  };
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
+  Uint FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::counter{0};
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
+  FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::FieldCollectionBase()
+    :id(counter++){}
+
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
+  void FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::register_field(Field_p &&field) {
+    auto&& search_it = this->fields.find(field->get_name());
+    auto&& does_exist = search_it != this->fields.end();
+    if (does_exist) {
+      std::stringstream err_str;
+      err_str << "a field named " << field->get_name()
+              << "is already registered in this field collection. "
+              << "Currently registered fields: ";
+      for (const auto& name_field_pair: this->fields) {
+        err_str << ", " << name_field_pair.first;
+      }
+      throw FieldCollectionError(err_str.str());
+    }
+    this->fields[field->get_name()] = std::move(field);
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
+  constexpr Dim_t FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
+  spatial_dim() {
+    return DimS;
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
+  Dim_t FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
+  get_spatial_dim() const {
+    return DimS;
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
+  typename FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::Field&
+  FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
+  operator[](std::string unique_name) {
+    return *(this->fields[unique_name]);
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
+  typename FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::Field&
+  FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
+  at(std::string unique_name) {
+    return *(this->fields.at(unique_name));
+  }
+
+
+}  // muSpectre
+
+#endif /* FIELD_COLLECTION_BASE_H */
diff --git a/src/common/field_collection_global.hh b/src/common/field_collection_global.hh
new file mode 100644
index 0000000..8289924
--- /dev/null
+++ b/src/common/field_collection_global.hh
@@ -0,0 +1,157 @@
+/**
+ * file   field_collection_global.hh
+ *
+ * @author Till Junge <till.junge@altermail.ch>
+ *
+ * @date   05 Nov 2017
+ *
+ * @brief  FieldCollection base-class for global fields
+ *
+ * @section LICENCE
+ *
+ * Copyright (C) 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.
+ */
+
+
+#ifndef FIELD_COLLECTION_GLOBAL_H
+#define FIELD_COLLECTION_GLOBAL_H
+
+#include "common/field_collection_base.hh"
+#include "common/ccoord_operations.hh"
+
+
+namespace muSpectre {
+
+  /* ---------------------------------------------------------------------- */
+  //! DimS spatial dimension (dimension of problem)
+  //! DimM material_dimension (dimension of constitutive law)
+  template<Dim_t DimS, Dim_t DimM>
+  class GlobalFieldCollection:
+    public FieldCollectionBase<DimS, DimM, GlobalFieldCollection<DimS, DimM>>
+  {
+  public:
+    using Parent = FieldCollectionBase
+      <DimS, DimM, GlobalFieldCollection<DimS, DimM>>;
+    using Ccoord = typename Parent::Ccoord;
+    using Field_p = typename Parent::Field_p;
+    //! Default constructor
+    GlobalFieldCollection();
+
+    //! Copy constructor
+    GlobalFieldCollection(const GlobalFieldCollection &other) = delete;
+
+    //! Move constructor
+    GlobalFieldCollection
+      (GlobalFieldCollection &&other) noexcept = delete;
+
+    //! Destructor
+    virtual ~GlobalFieldCollection() noexcept = default;
+
+    //! Copy assignment operator
+    GlobalFieldCollection&
+    operator=(const GlobalFieldCollection &other) = delete;
+
+    //! Move assignment operator
+    GlobalFieldCollection&
+    operator=(GlobalFieldCollection &&other) noexcept = delete;
+
+    /** allocate memory, etc. At this point, the collection is
+        informed aboud the size and shape of the domain (through the
+        sizes parameter). The job of initialise is to make sure that
+        all fields are either of size 0, in which case they need to be
+        allocated, or are of the same size as the product of 'sizes'
+        any field of a different size is wrong TODO: check whether it
+        makes sense to put a runtime check here
+     **/
+    inline void initialise(Ccoord sizes);
+
+    //! return the pixel sizes
+    inline const Ccoord & get_sizes() const;
+
+    //! returns the linear index corresponding to cell coordinates
+    inline size_t get_index(Ccoord && ccoord) const;
+    //! returns the cell coordinates corresponding to a linear index
+    inline Ccoord get_ccoord(size_t index) const;
+
+    static constexpr inline Dim_t spatial_dim() {return DimS;}
+    static constexpr inline Dim_t material_dim() {return DimM;}
+  protected:
+    //! number of discretisation cells in each of the DimS spatial directions
+    Ccoord sizes{0};
+  private:
+  };
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  GlobalFieldCollection<DimS, DimM>::GlobalFieldCollection()
+    :Parent(){}
+
+
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  void GlobalFieldCollection<DimS, DimM>::
+  initialise(Ccoord sizes) {
+    this->size_ = std::accumulate(sizes.begin(), sizes.end(), 1,
+                                   std::multiplies<Dim_t>());
+    this->sizes = sizes;
+    std::for_each(std::begin(this->fields), std::end(this->fields),
+                  [this](auto && item) {
+                    auto && field = *item.second;
+                    const auto field_size = field.size();
+                    if (field_size == 0) {
+                      field.resize(this->size());
+                    } else if (field_size != this->size()) {
+                      std::stringstream err_stream;
+                      err_stream << "Field '" << field.get_name()
+                                 << "' contains " << field_size
+                                 << " entries, but the field collection "
+                                 << "has " << this->size() << " pixels";
+                      throw FieldCollectionError(err_stream.str());
+                    }
+                  });
+    this->is_initialised = true;
+  }
+
+  //----------------------------------------------------------------------------//
+  //! return the pixel sizes
+  template <Dim_t DimS, Dim_t DimM>
+  const typename GlobalFieldCollection<DimS, DimM>::Ccoord &
+  GlobalFieldCollection<DimS, DimM>::get_sizes() const {
+    return this->sizes;
+  }
+
+  //----------------------------------------------------------------------------//
+  //! returns the cell coordinates corresponding to a linear index
+  template <Dim_t DimS, Dim_t DimM>
+  typename GlobalFieldCollection<DimS, DimM>::Ccoord
+  GlobalFieldCollection<DimS, DimM>::get_ccoord(size_t index) const {
+    return CcoordOps::get_ccoord(this->get_sizes(), std::move(index));
+  }
+
+  //----------------------------------------------------------------------------//
+  //! returns the linear index corresponding to cell coordinates
+  template <Dim_t DimS, Dim_t DimM>
+  size_t
+  GlobalFieldCollection<DimS, DimM>::get_index(Ccoord && ccoord) const {
+    return CcoordOps::get_index(this->get_sizes(), std::move(ccoord));
+  }
+
+}  // muSpectre
+
+#endif /* FIELD_COLLECTION_GLOBAL_H */
diff --git a/src/common/field_collection_local.hh b/src/common/field_collection_local.hh
new file mode 100644
index 0000000..7a84f40
--- /dev/null
+++ b/src/common/field_collection_local.hh
@@ -0,0 +1,160 @@
+/**
+ * file   field_collection_local.hh
+ *
+ * @author Till Junge <till.junge@altermail.ch>
+ *
+ * @date   05 Nov 2017
+ *
+ * @brief  FieldCollection base-class for local fields
+ *
+ * @section LICENCE
+ *
+ * Copyright (C) 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.
+ */
+
+#ifndef FIELD_COLLECTION_LOCAL_H
+#define FIELD_COLLECTION_LOCAL_H
+
+#include "common/field_collection_base.hh"
+
+namespace muSpectre {
+
+  /* ---------------------------------------------------------------------- */
+  //! DimS spatial dimension (dimension of problem)
+  //! DimM material_dimension (dimension of constitutive law)
+  template<Dim_t DimS, Dim_t DimM>
+  class LocalFieldCollection:
+    public FieldCollectionBase<DimS, DimM, LocalFieldCollection<DimS, DimM>>
+  {
+  public:
+    using Parent = FieldCollectionBase<DimS, DimM,
+                                       LocalFieldCollection<DimS, DimM>>;
+    using Ccoord = typename Parent::Ccoord;
+    using Field_p = typename Parent::Field_p;
+
+    //! Default constructor
+    LocalFieldCollection();
+
+    //! Copy constructor
+    LocalFieldCollection(const LocalFieldCollection &other) = delete;
+
+    //! Move constructor
+    LocalFieldCollection(LocalFieldCollection &&other) noexcept = delete;
+
+    //! Destructor
+    virtual ~LocalFieldCollection() noexcept = default;
+
+    //! Copy assignment operator
+    LocalFieldCollection& operator=(const LocalFieldCollection &other) = delete;
+
+    //! Move assignment operator
+    LocalFieldCollection& operator=(LocalFieldCollection &&other) noexcept = delete;
+
+    //! add a pixel/voxel to the field collection
+    inline void add_pixel(const Ccoord & local_ccoord);
+
+    /** allocate memory, etc. at this point, the field collection
+        knows how many entries it should have from the size of the
+        coords containes (which grows by one every time add_pixel is
+        called. The job of initialise is to make sure that all fields
+        are either of size 0, in which case they need to be allocated,
+        or are of the same size as the product of 'sizes' any field of
+        a different size is wrong TODO: check whether it makes sense
+        to put a runtime check here
+     **/
+    inline void initialise();
+
+
+    //! returns the linear index corresponding to cell coordinates
+    inline size_t get_index(Ccoord && ccoord) const;
+    //! returns the cell coordinates corresponding to a linear index
+    inline Ccoord get_ccoord(size_t index) const;
+
+
+  protected:
+    //! container of pixel coords for non-global collections
+    std::vector<Ccoord> ccoords;
+    //! container of indices for non-global collections (slow!)
+    std::map<Ccoord, size_t> indices;
+  private:
+  };
+
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  LocalFieldCollection<DimS, DimM>::LocalFieldCollection()
+    :Parent(){}
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  void LocalFieldCollection<DimS, DimM>::
+  add_pixel(const Ccoord & local_ccoord) {
+    if (this->is_initialised) {
+      throw FieldCollectionError
+        ("once a field collection has been initialised, you can't add new "
+         "pixels.");
+    }
+    this->indices[local_ccoord] = this->ccoords.size();
+    this->ccoords.push_back(local_ccoord);
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  void LocalFieldCollection<DimS, DimM>::
+  initialise() {
+    this->size_ = this->ccoords.size();
+    std::for_each(std::begin(this->fields), std::end(this->fields),
+                  [this](auto && item) {
+                    auto && field = *item.second;
+                    const auto field_size = field.size();
+                    if (field_size == 0) {
+                      field.resize(this->size());
+                    } else if (field_size != this->size()) {
+                      std::stringstream err_stream;
+                      err_stream << "Field '" << field.get_name()
+                                 << "' contains " << field_size
+                                 << " entries, but the field collection "
+                                 << "has " << this->size() << " pixels";
+                      throw FieldCollectionError(err_stream.str());
+                    }
+                  });
+    this->is_initialised = true;
+  }
+
+
+  //----------------------------------------------------------------------------//
+  //! returns the linear index corresponding to cell coordinates
+  template <Dim_t DimS, Dim_t DimM>
+  size_t
+  LocalFieldCollection<DimS, DimM>::get_index(Ccoord && ccoord) const {
+    return this->indices[std::move(ccoord)];
+  }
+
+
+  //----------------------------------------------------------------------------//
+  //! returns the cell coordinates corresponding to a linear index
+  template <Dim_t DimS, Dim_t DimM>
+  typename LocalFieldCollection<DimS, DimM>::Ccoord
+  LocalFieldCollection<DimS, DimM>::get_ccoord(size_t index) const {
+    return this->ccoords[std::move(index)];
+  }
+
+
+}  // muSpectre
+
+#endif /* FIELD_COLLECTION_LOCAL_H */
diff --git a/src/materials/material.cc b/src/materials/material.cc
index dccd629..a228564 100644
--- a/src/materials/material.cc
+++ b/src/materials/material.cc
@@ -1,46 +1,46 @@
 /**
  * file   material.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   01 Nov 2017
  *
  * @brief  implementation of materi
  *
  * @section LICENCE
  *
  * Copyright (C) 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 "material.hh"
+#include "materials/material.hh"
 
 namespace muSpectre {
 
   //----------------------------------------------------------------------------//
   template<Dim_t s_dim, Dim_t m_dim>
   MaterialBase<s_dim, m_dim>::MaterialBase(std::string name):name(name){
     static_assert((m_dim == oneD)||
                   (m_dim == twoD)||
                   (m_dim == threeD), "only 1, 2, or threeD supported");
   }
 
   template class MaterialBase<2, 2>;
   template class MaterialBase<2, 3>;
   template class MaterialBase<3, 3>;
 
 }  // muSpectre
diff --git a/tests/test_field_collections.cc b/tests/test_field_collections.cc
index 43fdd61..09f75b2 100644
--- a/tests/test_field_collections.cc
+++ b/tests/test_field_collections.cc
@@ -1,547 +1,548 @@
 /**
  * file   test_field_collections.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   20 Sep 2017
  *
  * @brief  Test the FieldCollection classes which provide fast optimized iterators
  *         over run-time typed fields
  *
  * @section LICENCE
  *
  * Copyright (C) 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 <stdexcept>
 #include <boost/mpl/list.hpp>
 #include <random>
 #include <type_traits>
 #include <sstream>
 #include <string>
 
 #include "common/common.hh"
 #include "common/ccoord_operations.hh"
 #include "common/test_goodies.hh"
 #include "tests.hh"
 #include "common/field_collection.hh"
 #include "common/field.hh"
 #include "common/field_map.hh"
 
 namespace muSpectre {
 
   BOOST_AUTO_TEST_SUITE(field_collection_tests);
 
   //! Test fixture for simple tests on single field in collection
   template <Dim_t DimS, Dim_t DimM, bool Global>
   struct FC_fixture:
     public FieldCollection<DimS, DimM, Global> {
     FC_fixture()
       :fc() {}
     inline static constexpr Dim_t sdim(){return DimS;}
     inline static constexpr Dim_t mdim(){return DimM;}
     inline static constexpr bool global(){return Global;}
     using FC_t = FieldCollection<DimS, DimM, Global>;
     FC_t fc;
   };
 
   using test_collections = boost::mpl::list<FC_fixture<2, 2, true>,
                                             FC_fixture<2, 3, true>,
                                             FC_fixture<3, 3, true>,
                                             FC_fixture<2, 2, false>,
                                             FC_fixture<2, 3, false>,
                                             FC_fixture<3, 3, false>>;
 
   BOOST_AUTO_TEST_CASE(simple) {
     const Dim_t sdim = 2;
     const Dim_t mdim = 2;
     using FC_t = FieldCollection<sdim, mdim>;
     FC_t fc;
 
     BOOST_CHECK_EQUAL(FC_t::spatial_dim(), sdim);
     BOOST_CHECK_EQUAL(fc.get_spatial_dim(), sdim);
   }
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(Simple_construction_test, F, test_collections, F) {
     BOOST_CHECK_EQUAL(F::FC_t::spatial_dim(), F::sdim());
     BOOST_CHECK_EQUAL(F::fc.get_spatial_dim(), F::sdim());
 
   }
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(get_field2_test, F, test_collections, F) {
     const auto order{2};
 
     using FC_t = typename F::FC_t;
     using TF_t = TensorField<FC_t, Real, order, F::mdim()>;
     auto && myfield = make_field<TF_t>("TensorField real 2", F::fc);
 
     using TensorMap = TensorFieldMap<FC_t, Real, order,F::mdim()>;
     using MatrixMap = MatrixFieldMap<FC_t, Real, F::mdim(), F::mdim()>;
     using ArrayMap = ArrayFieldMap<FC_t, Real, F::mdim(), F::mdim()>;
 
     TensorMap TFM(myfield);
     MatrixMap MFM(myfield);
     ArrayMap  AFM(myfield);
 
     BOOST_CHECK_EQUAL(TFM.info_string(),
                       "Tensor(d, "+ std::to_string(order) +
                       "_o, " + std::to_string(F::mdim()) + "_d)");
     BOOST_CHECK_EQUAL(MFM.info_string(),
                       "Matrix(d, "+ std::to_string(F::mdim()) +
                       "x" + std::to_string(F::mdim()) + ")");
     BOOST_CHECK_EQUAL(AFM.info_string(),
                       "Array(d, "+ std::to_string(F::mdim()) +
                       "x" + std::to_string(F::mdim()) + ")");
 
   }
 
   constexpr Dim_t order{4}, matrix_order{2};
   //! Test fixture for multiple fields in the collection
   template <Dim_t DimS, Dim_t DimM, bool Global>
   struct FC_multi_fixture{
     FC_multi_fixture()
       :fc() {
       //add Real tensor field
       make_field<TensorField<FC_t, Real, order, DimM>>
         ("Tensorfield Real o4", fc);
       //add Real tensor field
       make_field<TensorField<FC_t, Real, matrix_order, DimM>>
         ("Tensorfield Real o2", fc);
       //add integer scalar field
       make_field<ScalarField<FC_t, Int>>
         ("integer Scalar", fc);
       //add complex matrix field
       make_field<MatrixField<FC_t, Complex, DimS, DimM>>
         ("Matrixfield Complex sdim x mdim", fc);
 
     }
     inline static constexpr Dim_t sdim(){return DimS;}
     inline static constexpr Dim_t mdim(){return DimM;}
     inline static constexpr bool global(){return Global;}
     using FC_t = FieldCollection<DimS, DimM, Global>;
     FC_t fc;
   };
 
   using mult_collections = boost::mpl::list<FC_multi_fixture<2, 2, true>,
                                             FC_multi_fixture<2, 3, true>,
                                             FC_multi_fixture<3, 3, true>,
                                             FC_multi_fixture<2, 2, false>,
                                             FC_multi_fixture<2, 3, false>,
                                             FC_multi_fixture<3, 3, false>>;
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(multi_field_test, F, mult_collections, F) {
     using FC_t = typename F::FC_t;
     // possible maptypes for Real tensor fields
     using T_type = Real;
     using T_TFM1_t = TensorFieldMap<FC_t, T_type, order, F::mdim()>;
     using T_TFM2_t = TensorFieldMap<FC_t, T_type, 2, F::mdim()*F::mdim()>; //! dangerous
     // impossible maptypes for Real tensor fields
     using T_SFM_t = ScalarFieldMap<FC_t, T_type>;
     using T_MFM_t = MatrixFieldMap<FC_t, T_type, 1, 1>;
     using T_AFM_t = ArrayFieldMap<FC_t, T_type, 1, 1>;
     using T_MFMw1_t = MatrixFieldMap<FC_t, Int, 1, 2>;
     using T_MFMw2_t = MatrixFieldMap<FC_t, Real, 1, 2>;
     using T_MFMw3_t = MatrixFieldMap<FC_t, Complex, 1, 2>;
     const std::string T_name{"Tensorfield Real o4"};
     const std::string T_name_w{"TensorField Real o4 wrongname"};
 
     BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name)), FieldInterpretationError);
     BOOST_CHECK_NO_THROW(T_TFM1_t(F::fc.at(T_name)));
     BOOST_CHECK_NO_THROW(T_TFM2_t(F::fc.at(T_name)));
     BOOST_CHECK_THROW(T_MFM_t(F::fc.at(T_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(T_AFM_t(F::fc.at(T_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(T_MFMw1_t(F::fc.at(T_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(T_MFMw3_t(F::fc.at(T_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name_w)), std::out_of_range);
 
     // possible maptypes for integer scalar fields
     using S_type = Int;
     using S_SFM_t = ScalarFieldMap<FC_t, S_type>;
     using S_TFM1_t = TensorFieldMap<FC_t, S_type, 1, 1>;
     using S_TFM2_t = TensorFieldMap<FC_t, S_type, 2, 1>;
     using S_MFM_t = MatrixFieldMap<FC_t, S_type, 1, 1>;
     using S_AFM_t = ArrayFieldMap<FC_t, S_type, 1, 1>;
     // impossible maptypes for integer scalar fields
     using S_MFMw1_t = MatrixFieldMap<FC_t, Int, 1, 2>;
     using S_MFMw2_t = MatrixFieldMap<FC_t, Real, 1, 2>;
     using S_MFMw3_t = MatrixFieldMap<FC_t, Complex, 1, 2>;
     const std::string S_name{"integer Scalar"};
     const std::string S_name_w{"integer Scalar wrongname"};
 
     BOOST_CHECK_NO_THROW(S_SFM_t(F::fc.at(S_name)));
     BOOST_CHECK_NO_THROW(S_TFM1_t(F::fc.at(S_name)));
     BOOST_CHECK_NO_THROW(S_TFM2_t(F::fc.at(S_name)));
     BOOST_CHECK_NO_THROW(S_MFM_t(F::fc.at(S_name)));
     BOOST_CHECK_NO_THROW(S_AFM_t(F::fc.at(S_name)));
     BOOST_CHECK_THROW(S_MFMw1_t(F::fc.at(S_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(S_MFMw3_t(F::fc.at(S_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(S_SFM_t(F::fc.at(S_name_w)), std::out_of_range);
 
     // possible maptypes for complex matrix fields
     using M_type = Complex;
     using M_MFM_t = MatrixFieldMap<FC_t, M_type, F::sdim(), F::mdim()>;
     using M_AFM_t = ArrayFieldMap<FC_t, M_type, F::sdim(), F::mdim()>;
     // impossible maptypes for complex matrix fields
     using M_SFM_t = ScalarFieldMap<FC_t, M_type>;
     using M_MFMw1_t = MatrixFieldMap<FC_t, Int, 1, 2>;
     using M_MFMw2_t = MatrixFieldMap<FC_t, Real, 1, 2>;
     using M_MFMw3_t = MatrixFieldMap<FC_t, Complex, 1, 2>;
     const std::string M_name{"Matrixfield Complex sdim x mdim"};
     const std::string M_name_w{"Matrixfield Complex sdim x mdim wrongname"};
 
     BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name)), FieldInterpretationError);
     BOOST_CHECK_NO_THROW(M_MFM_t(F::fc.at(M_name)));
     BOOST_CHECK_NO_THROW(M_AFM_t(F::fc.at(M_name)));
     BOOST_CHECK_THROW(M_MFMw1_t(F::fc.at(M_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(M_MFMw3_t(F::fc.at(M_name)), FieldInterpretationError);
     BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name_w)), std::out_of_range);
 
   }
 
   /* ---------------------------------------------------------------------- */
   //! Check whether fields can be initialized
   using mult_collections_t = boost::mpl::list<FC_multi_fixture<2, 2, true>,
                                               FC_multi_fixture<2, 3, true>,
                                               FC_multi_fixture<3, 3, true>>;
   using mult_collections_f = boost::mpl::list<FC_multi_fixture<2, 2, false>,
                                               FC_multi_fixture<2, 3, false>,
                                               FC_multi_fixture<3, 3, false>>;
 
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_glob, F, mult_collections_t, F) {
     Ccoord_t<F::sdim()> size;
     for (auto && s: size) {
       s = 3;
     }
     BOOST_CHECK_NO_THROW(F::fc.initialise(size));
   }
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca, F, mult_collections_f, F) {
     testGoodies::RandRange<Int> rng;
     for (int i = 0; i < 7; ++i) {
       Ccoord_t<F::sdim()> pixel;
       for (auto && s: pixel) {
         s = rng.randval(0, 7);
       }
       F::fc.add_pixel(pixel);
     }
 
     BOOST_CHECK_NO_THROW(F::fc.initialise());
   }
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca_with_push_back, F, mult_collections_f, F) {
     constexpr auto mdim = F::mdim();
     constexpr int nb_pix = 7;
     testGoodies::RandRange<Int> rng;
     using ftype = internal::TypedFieldBase<decltype(F::fc), Real,  mdim*mdim*mdim*mdim>;
     using stype = Eigen::Array<Real, mdim*mdim*mdim*mdim, 1>;
     auto & field = reinterpret_cast<ftype&>(F::fc["Tensorfield Real o4"]);
     field.push_back(stype());
     for (int i = 0; i < nb_pix; ++i) {
       Ccoord_t<F::sdim()> pixel;
       for (auto && s: pixel) {
         s = rng.randval(0, 7);
       }
       F::fc.add_pixel(pixel);
     }
 
     BOOST_CHECK_THROW(F::fc.initialise(), FieldCollectionError);
     for (int i = 0; i < nb_pix-1; ++i) {
       field.push_back(stype());
     }
     BOOST_CHECK_NO_THROW(F::fc.initialise());
 
   }
 
   //! Test fixture for iterators over multiple fields
   template <Dim_t DimS, Dim_t DimM, bool Global>
   struct FC_iterator_fixture
     : public FC_multi_fixture<DimS, DimM, Global> {
     using Parent = FC_multi_fixture<DimS, DimM, Global>;
     FC_iterator_fixture()
       :Parent() {
       this-> fill();
     }
 
     template <bool isGlobal = Global>
     std::enable_if_t<isGlobal> fill() {
       static_assert(Global==isGlobal, "You're breaking my SFINAE plan");
       Ccoord_t<Parent::sdim()> size;
       for (auto && s: size) {
         s = cube_size();
       }
       this->fc.initialise(size);
     }
 
     template <bool notGlobal = !Global>
     std::enable_if_t<notGlobal> fill (int dummy = 0) {
       static_assert(notGlobal != Global, "You're breaking my SFINAE plan");
       testGoodies::RandRange<Int> rng;
+      this->fc.add_pixel({0,0});
       for (int i = 0*dummy; i < sele_size(); ++i) {
         Ccoord_t<Parent::sdim()> pixel;
         for (auto && s: pixel) {
           s = rng.randval(0, 7);
         }
         this->fc.add_pixel(pixel);
       }
 
       this->fc.initialise();
     }
 
     constexpr static Dim_t cube_size() {return 3;}
     constexpr static Dim_t sele_size() {return 7;}
   };
 
   using iter_collections = boost::mpl::list<FC_iterator_fixture<2, 2, true>,
                                             FC_iterator_fixture<2, 3, true>,
                                             FC_iterator_fixture<3, 3, true>,
                                             FC_iterator_fixture<2, 2, false>,
                                             FC_iterator_fixture<2, 3, false>,
                                             FC_iterator_fixture<3, 3, false>>;
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(iter_field_test, F, iter_collections, F) {
     using FC_t = typename F::Parent::FC_t;
     using Tensor4Map = TensorFieldMap<FC_t, Real, order, F::Parent::mdim()>;
     Tensor4Map T4map{F::fc["Tensorfield Real o4"]};
     F::fc["Tensorfield Real o4"].set_zero();
     for (auto && tens:T4map) {
       BOOST_CHECK_EQUAL(Real(Eigen::Tensor<Real, 0>(tens.abs().sum().eval())()), 0);
     }
 
     using Tensor2Map = TensorFieldMap<FC_t, Real, matrix_order, F::Parent::mdim()>;
     using MSqMap = MatrixFieldMap<FC_t, Real, F::Parent::mdim(), F::Parent::mdim()>;
     using ASqMap =  ArrayFieldMap<FC_t, Real, F::Parent::mdim(), F::Parent::mdim()>;
     Tensor2Map T2map{F::fc["Tensorfield Real o2"]};
     MSqMap Mmap{F::fc["Tensorfield Real o2"]};
     ASqMap Amap{F::fc["Tensorfield Real o2"]};
     auto t2_it = T2map.begin();
     auto t2_it_end = T2map.end();
     auto m_it = Mmap.begin();
     auto a_it = Amap.begin();
     for (; t2_it != t2_it_end; ++t2_it, ++m_it, ++a_it) {
       t2_it->setRandom();
       auto && m = *m_it;
       bool comp = (m == a_it->matrix());
       BOOST_CHECK(comp);
     }
 
     using ScalarMap = ScalarFieldMap<FC_t, Int>;
     ScalarMap s_map{F::fc["integer Scalar"]};
     for (Uint i = 0; i < s_map.size(); ++i) {
       s_map[i] = i;
     }
     Uint counter{0};
 
     for (const auto& val: s_map) {
       BOOST_CHECK_EQUAL(counter++, val);
     }
 
   }
 
   using glob_iter_colls = boost::mpl::list<FC_iterator_fixture<2, 2, true>,
                                            FC_iterator_fixture<2, 3, true>,
                                            FC_iterator_fixture<3, 3, true>>;
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(ccoord_indexing_test, F, glob_iter_colls, F) {
     using FC_t = typename F::Parent::FC_t;
     using ScalarMap = ScalarFieldMap<FC_t, Int>;
     ScalarMap s_map{F::fc["integer Scalar"]};
     for (Uint i = 0; i < s_map.size(); ++i) {
       s_map[i] = i;
     }
 
 
     for (size_t i = 0; i < CcoordOps::get_size(F::fc.get_sizes()); ++i) {
       Ccoord_t<F::sdim()> sizes(F::fc.get_sizes());
       Ccoord_t<F::sdim()> ccoords(CcoordOps::get_ccoord(sizes, i));
       BOOST_CHECK_EQUAL(CcoordOps::get_index(F::fc.get_sizes(),
                                              CcoordOps::get_ccoord(F::fc.get_sizes(), i)), i);
     }
 
   }
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(iterator_methods_test, F, iter_collections, F) {
     using FC_t = typename F::Parent::FC_t;
     using Tensor4Map = TensorFieldMap<FC_t, Real, order, F::Parent::mdim()>;
     Tensor4Map T4map{F::fc["Tensorfield Real o4"]};
     using it_t = typename Tensor4Map::iterator;
     std::ptrdiff_t diff{3}; // arbitrary, as long as it is smaller than the container size
 
     // check constructors
     auto itstart = T4map.begin(); // standard way of obtaining iterator
     auto itend = T4map.end(); // ditto
 
     it_t it1{T4map};
     it_t it2{T4map, false};
     it_t it3{T4map, size_t(diff)};
     BOOST_CHECK(itstart == itstart);
     BOOST_CHECK(itstart != itend);
     BOOST_CHECK_EQUAL(itstart, it1);
     BOOST_CHECK_EQUAL(itend, it2);
 
     // check ostream operator
     std::stringstream response;
     response << it3;
     BOOST_CHECK_EQUAL
       (response.str(),
        std::string
        ("iterator on field 'Tensorfield Real o4', entry ") +
        std::to_string(diff));
 
     // check copy, move, and assigment constructor (and operator+)
     it_t itcopy = it1;
     it_t itmove = std::move(T4map.begin());
     it_t it4 = it1+diff;
     it_t it5(it2);
     it_t tmp(it2);
     it_t it6(std::move(tmp));
     it_t it7 = it4 -diff;
     BOOST_CHECK_EQUAL(itcopy, it1);
     BOOST_CHECK_EQUAL(itmove, it1);
     BOOST_CHECK_EQUAL(it4, it3);
     BOOST_CHECK_EQUAL(it5, it2);
     BOOST_CHECK_EQUAL(it6, it5);
     BOOST_CHECK_EQUAL(it7, it1);
 
     // check increments/decrements
     BOOST_CHECK_EQUAL(it1++, itcopy);    // post-increment
     BOOST_CHECK_EQUAL(it1, itcopy+1);
     BOOST_CHECK_EQUAL(--it1, itcopy);    // pre -decrement
     BOOST_CHECK_EQUAL(++it1, itcopy+1);  // pre -increment
     BOOST_CHECK_EQUAL(it1--, itcopy+1);  // post-decrement
     BOOST_CHECK_EQUAL(it1, itcopy);
 
 
     // dereference and member-of-pointer check
     Eigen::Tensor<Real, 4> Tens = *it1;
     Eigen::Tensor<Real, 4> Tens2 = *itstart;
     Eigen::Tensor<bool, 0> check = (Tens==Tens2).all();
     BOOST_CHECK_EQUAL(bool(check()), true);
 
     BOOST_CHECK_NO_THROW(itstart->setZero());
 
     //check access subscripting
     auto T3a = *it3;
     auto T3b = itstart[diff];
     BOOST_CHECK(bool(Eigen::Tensor<bool, 0>((T3a==T3b).all())()));
 
     // div. comparisons
     BOOST_CHECK_LT(itstart, itend);
     BOOST_CHECK(!(itend < itstart));
     BOOST_CHECK(!(itstart < itstart));
 
     BOOST_CHECK_LE(itstart, itend);
     BOOST_CHECK_LE(itstart, itstart);
     BOOST_CHECK(!(itend <= itstart));
 
     BOOST_CHECK_GT(itend, itstart);
     BOOST_CHECK(!(itend>itend));
     BOOST_CHECK(!(itstart>itend));
 
     BOOST_CHECK_GE(itend, itstart);
     BOOST_CHECK_GE(itend, itend);
     BOOST_CHECK(!(itstart >= itend));
 
     // check assignment increment/decrement
     BOOST_CHECK_EQUAL(it1+=diff, it3);
     BOOST_CHECK_EQUAL(it1-=diff, itstart);
 
     // check cell coordinates
     using Ccoord = Ccoord_t<F::sdim()>;
     Ccoord a{itstart.get_ccoord()};
-    Ccoord b{};
+    Ccoord b{0};
 
-    // Weirdly, boost::has_left_shift<std::ostream, T>::value is falso for Ccoords, even though the operator is implemented :(
+    // Weirdly, boost::has_left_shift<std::ostream, T>::value is false for Ccoords, even though the operator is implemented :(
     //BOOST_CHECK_EQUAL(a, b);
     bool check2 = (a==b);
     BOOST_CHECK(check2);
   }
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_tensor_iter_test, F, iter_collections, F) {
     using FC_t = typename F::Parent::FC_t;
     using Tensor4Map = TensorFieldMap<FC_t, Real, order, F::Parent::mdim()>;
     Tensor4Map T4map{F::fc["Tensorfield Real o4"]};
 
     using T_t = typename Tensor4Map::T_t;
     Eigen::TensorMap<const T_t> Tens2(T4map[0].data(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim());
 
     for (auto it = T4map.cbegin(); it != T4map.cend(); ++it) {
       // maps to const tensors can't be initialised with a const pointer this sucks
       auto&&  tens = *it;
       auto&&  ptr = tens.data();
 
       static_assert(std::is_pointer<std::remove_reference_t<decltype(ptr)>>::value, "should be getting a pointer");
       //static_assert(std::is_const<std::remove_pointer_t<decltype(ptr)>>::value, "should be const");
 
       // If Tensor were written well, above static_assert should pass, and the
       // following check shouldn't. If it get triggered, it means that a newer
       // version of Eigen now does have const-correct
       // TensorMap<const Tensor<...>. This means that const-correct field maps
       // are then also possible for tensors
       BOOST_CHECK(!std::is_const<std::remove_pointer_t<decltype(ptr)>>::value);
 
     }
   }
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_matrix_iter_test, F, iter_collections, F) {
     using FC_t = typename F::Parent::FC_t;
     using MatrixMap = MatrixFieldMap<FC_t, Complex, F::sdim(), F::mdim()>;
     MatrixMap Mmap{F::fc["Matrixfield Complex sdim x mdim"]};
 
     for (auto it = Mmap.cbegin(); it != Mmap.cend(); ++it) {
       // maps to const tensors can't be initialised with a const pointer this sucks
       auto&&  mat = *it;
       auto&&  ptr = mat.data();
 
       static_assert(std::is_pointer<std::remove_reference_t<decltype(ptr)>>::value, "should be getting a pointer");
       //static_assert(std::is_const<std::remove_pointer_t<decltype(ptr)>>::value, "should be const");
 
       // If Matrix were written well, above static_assert should pass, and the
       // following check shouldn't. If it get triggered, it means that a newer
       // version of Eigen now does have const-correct
       // MatrixMap<const Matrix<...>. This means that const-correct field maps
       // are then also possible for matrices
       BOOST_CHECK(!std::is_const<std::remove_pointer_t<decltype(ptr)>>::value);
 
     }
   }
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_scalar_iter_test, F, iter_collections, F) {
     using FC_t = typename F::Parent::FC_t;
     using ScalarMap = ScalarFieldMap<FC_t, Int>;
     ScalarMap Smap{F::fc["integer Scalar"]};
 
     for (auto it = Smap.cbegin(); it != Smap.cend(); ++it) {
       auto&& scal = *it;
       static_assert(std::is_const<std::remove_reference_t<decltype(scal)>>::value,
                     "referred type should be const");
       static_assert(std::is_lvalue_reference<decltype(scal)>::value,
                     "Should have returned an lvalue ref");
 
     }
   }
 
   BOOST_AUTO_TEST_SUITE_END();
 
 }  // muSpectre