diff --git a/src/common/field.hh b/src/common/field.hh
index a6802f3..dda09be 100644
--- a/src/common/field.hh
+++ b/src/common/field.hh
@@ -1,535 +1,536 @@
 /**
  * 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 <sstream>
 #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:
       using collection_t = FieldCollection;
 
       //! 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;
+      friend typename FieldCollection::Parent;
 
     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, bool isConst>
     class FieldMap;
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents,
               bool ArrayStore=false>
     class TypedFieldBase: public FieldBase<FieldCollection>
     {
       friend class FieldMap<FieldCollection, T, NbComponents, true>;
       friend class FieldMap<FieldCollection, T, NbComponents, false>;
     public:
       constexpr static auto nb_components{NbComponents};
       using Parent = FieldBase<FieldCollection>;
       using Parent::collection_t;
       using Scalar = T;
       using Base = Parent;
       //using storage_type = Eigen::Array<T, Eigen::Dynamic, NbComponents>;
       using StoredType = Eigen::Array<T, NbComponents, 1>;
       using StorageType = std::conditional_t
         <ArrayStore,
          std::vector<StoredType,
                      Eigen::aligned_allocator<StoredType>>,
          std::vector<T,Eigen::aligned_allocator<T>>>;
       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;
       template <bool isArrayStore = ArrayStore>
       inline void push_back(const
                             std::enable_if_t<isArrayStore, StoredType> &
                             value);
       template <bool componentStore = !ArrayStore>
       inline std::enable_if_t<componentStore> push_back(const StoredType & value);
       size_t size() const override final;
 
       static TypedFieldBase & check_ref(Parent & other);
       static const TypedFieldBase & check_ref(const Base & parent);
 
       inline T* data() {return this->get_ptr_to_entry(0);}
       inline const T* data() const {return this->get_ptr_to_entry(0);}
     protected:
 
       template <bool isArray=ArrayStore>
       inline std::enable_if_t<isArray, T*> get_ptr_to_entry(const size_t&& index);
 
       template <bool isArray=ArrayStore>
       inline std::enable_if_t<isArray, const T*>
       get_ptr_to_entry(const size_t&& index) const;
 
       template <bool noArray = !ArrayStore>
       inline T* get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index);
 
       template <bool noArray = !ArrayStore>
       inline const T*
       get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) const;
 
       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 = 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 FieldType& make_field(std::string unique_name,
                              CollectionType & collection,
                              Args&&... args);
 
     static TensorField & check_ref(Base & other) {
       return static_cast<TensorField &>(Parent::check_ref(other));}
     static const TensorField & check_ref(const Base & other) {
       return static_cast<const TensorField &>(Parent::check_ref(other));}
 
   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 FieldType&  make_field(std::string unique_name,
                              CollectionType & collection,
                              Args&&... args);
 
     static MatrixField & check_ref(Base & other) {
       return static_cast<MatrixField &>(Parent::check_ref(other));}
     static const MatrixField & check_ref(const Base & other) {
       return static_cast<const MatrixField &>(Parent::check_ref(other));}
 
 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, bool ArrayStore>
     TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     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, bool ArrayStore>
     const std::type_info & TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     get_stored_typeid() const {
       return typeid(T);
     }
 
     /* ---------------------------------------------------------------------- */
     template<class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     void
     TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     set_zero() {
       std::fill(this->array.begin(), this->array.end(), T{});
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     size_t TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     size() const {
       if (ArrayStore) {
         return this->array.size();
       } else  {
         return this->array.size()/NbComponents;
       }
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
     TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::check_ref(Base & other) {
       if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) {
         std::string err ="Cannot create a Reference of requested type " +(
            "for field '" + other.get_name() + "' of type '" +
            other.get_stored_typeid().name() + "'");
         throw std::runtime_error
           (err);
       }
       //check size compatibility
       if (NbComponents != other.get_nb_components()) {
         throw std::runtime_error
           ("Cannot create a Reference to a field with " +
            std::to_string(NbComponents) + " components " +
            "for field '" + other.get_name() + "' with " +
            std::to_string(other.get_nb_components()) + " components");
       }
       return static_cast<TypedFieldBase&>(other);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     const TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
     TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     check_ref(const Base & other) {
       if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) {
         std::stringstream err_str{};
         err_str << "Cannot create a Reference of requested type "
                 << "for field '"  << other.get_name() << "' of type '"
                 << other.get_stored_typeid().name() << "'";
         throw std::runtime_error
           (err_str.str());
       }
       //check size compatibility
       if (NbComponents != other.get_nb_components()) {
         throw std::runtime_error
           ("Cannot create a Reference to a field with " +
            std::to_string(NbComponents) + " components " +
            "for field '" + other.get_name() + "' with " +
            std::to_string(other.get_nb_components()) + " components");
       }
       return static_cast<const TypedFieldBase&>(other);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool isArray>
     std::enable_if_t<isArray, T*>
     TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     get_ptr_to_entry(const size_t&& index) {
       static_assert (isArray == ArrayStore, "SFINAE");
       return &this->array[std::move(index)](0, 0);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool noArray>
     T* TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) {
       static_assert (noArray != ArrayStore, "SFINAE");
       return &this->array[NbComponents*std::move(index)];
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool isArray>
     std::enable_if_t<isArray, const T*>
     TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     get_ptr_to_entry(const size_t&& index) const {
       static_assert (isArray == ArrayStore, "SFINAE");
       return &this->array[std::move(index)](0, 0);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool noArray>
     const T* TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) const {
       static_assert (noArray != ArrayStore, "SFINAE");
       return &this->array[NbComponents*std::move(index)];
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     void TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     resize(size_t size) {
       if (ArrayStore) {
         this->array.resize(size);
       } else {
         this->array.resize(size*NbComponents);
       }
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool isArrayStore>
     void
     TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     push_back(const std::enable_if_t<isArrayStore,StoredType> & value) {
       static_assert(isArrayStore == ArrayStore, "SFINAE");
       this->array.push_back(value);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
     template <bool componentStore>
     std::enable_if_t<componentStore>
     TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
     push_back(const StoredType & value) {
       static_assert(componentStore != ArrayStore, "SFINAE");
       for (Dim_t i = 0; i < NbComponents; ++i) {
         this->array.push_back(value(i));
       }
     }
 
   }  // 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>
   FieldType &
   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_base.hh b/src/common/field_collection_base.hh
index 1c16da9..669c6cc 100644
--- a/src/common/field_collection_base.hh
+++ b/src/common/field_collection_base.hh
@@ -1,185 +1,188 @@
 /**
  * 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
     //! as 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());
     }
+    if (this->is_initialised) {
+      field->resize(this->size());
+    }
     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
index d04e0c5..75fd5db 100644
--- a/src/common/field_collection_global.hh
+++ b/src/common/field_collection_global.hh
@@ -1,186 +1,189 @@
 /**
  * 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;
     using iterator = typename CcoordOps::Pixels<DimS>::iterator;
     //! 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
     template <class CcoordRef>
     inline size_t get_index(CcoordRef && ccoord) const;
     //! returns the cell coordinates corresponding to a linear index
     inline Ccoord get_ccoord(size_t index) const;
 
     inline iterator begin();
     inline iterator end();
 
 
     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};
     CcoordOps::Pixels<DimS> pixels;
   private:
   };
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   GlobalFieldCollection<DimS, DimM>::GlobalFieldCollection()
     :Parent(), pixels{{0}}{}
 
 
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void GlobalFieldCollection<DimS, DimM>::
   initialise(Ccoord sizes) {
+    if (this->is_initialised) {
+      throw std::runtime_error("double initialisation");
+    }
     this->pixels = CcoordOps::Pixels<DimS>(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));
   }
 
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   typename GlobalFieldCollection<DimS, DimM>::iterator
   GlobalFieldCollection<DimS, DimM>::begin() {
     return this->pixels.begin();
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   typename GlobalFieldCollection<DimS, DimM>::iterator
   GlobalFieldCollection<DimS, DimM>::end() {
     return this->pixels.end();
   }
   //----------------------------------------------------------------------------//
   //! returns the linear index corresponding to cell coordinates
   template <Dim_t DimS, Dim_t DimM>
   template <class CcoordRef>
   size_t
   GlobalFieldCollection<DimS, DimM>::get_index(CcoordRef && ccoord) const {
     static_assert(std::is_same<
                     Ccoord,
                     std::remove_const_t<
                       std::remove_reference_t<CcoordRef>>>::value,
                   "can only be called with values or references of Ccoord");
     return CcoordOps::get_index(this->get_sizes(),
                                 std::forward<CcoordRef>(ccoord));
   }
 
 }  // muSpectre
 
 #endif /* FIELD_COLLECTION_GLOBAL_H */
diff --git a/src/common/field_collection_local.hh b/src/common/field_collection_local.hh
index 2348238..a8ff2ad 100644
--- a/src/common/field_collection_local.hh
+++ b/src/common/field_collection_local.hh
@@ -1,171 +1,174 @@
 /**
  * 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;
     using ccoords_container = std::vector<Ccoord>;
     using iterator = typename ccoords_container::iterator;
 
     //! 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
     template <class CcoordRef>
     inline size_t get_index(CcoordRef && ccoord) const;
     //! returns the cell coordinates corresponding to a linear index
     inline Ccoord get_ccoord(size_t index) const;
 
     inline iterator begin() {return this->ccoords.begin();}
     inline iterator end() {return this->ccoords.end();}
 
   protected:
     //! container of pixel coords for non-global collections
     ccoords_container 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() {
+    if (this->is_initialised) {
+      throw std::runtime_error("double initialisation");
+    }
     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>
   template <class CcoordRef>
   size_t
   LocalFieldCollection<DimS, DimM>::get_index(CcoordRef && ccoord) const {
     static_assert(std::is_same<
                     Ccoord,
                     std::remove_const_t<
                       std::remove_reference_t<CcoordRef>>>::value,
                   "can only be called with values or references of Ccoord");
     return this->indices.at(std::forward<CcoordRef>(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/fft/projection_base.hh b/src/fft/projection_base.hh
index 8fc8735..479372b 100644
--- a/src/fft/projection_base.hh
+++ b/src/fft/projection_base.hh
@@ -1,93 +1,100 @@
 /**
  * file   projection_base.hh
  *
  * @author Till Junge <till.junge@altermail.ch>
  *
  * @date   03 Dec 2017
  *
  * @brief  Base class for Projection operators
  *
  * @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 PROJECTION_BASE_H
 #define PROJECTION_BASE_H
 
 #include "common/common.hh"
 #include "common/field_collection.hh"
 #include "common/field.hh"
 #include "fft/fft_engine_base.hh"
 
 namespace muSpectre {
 
   template<class Projection>
   struct Projection_traits {
   };
 
   template <Dim_t DimS, Dim_t DimM>
   class ProjectionBase
   {
   public:
     using FFT_Engine = FFT_Engine_base<DimS, DimM>;
     using Ccoord = typename FFT_Engine::Ccoord;
+    using Rcoord = typename FFT_Engine::Rcoord;
     using GFieldCollection_t = typename FFT_Engine::GFieldCollection_t;
     using LFieldCollection_t = typename FFT_Engine::LFieldCollection_t;
     using Field_t = typename FFT_Engine::Field_t;
     using iterator = typename FFT_Engine::iterator;
 
     //! Default constructor
     ProjectionBase() = delete;
 
     //! Constructor with system sizes
     ProjectionBase(FFT_Engine & engine);
 
     //! Copy constructor
     ProjectionBase(const ProjectionBase &other) = delete;
 
     //! Move constructor
     ProjectionBase(ProjectionBase &&other) noexcept = default;
 
     //! Destructor
     virtual ~ProjectionBase() noexcept = default;
 
     //! Copy assignment operator
     ProjectionBase& operator=(const ProjectionBase &other) = delete;
 
     //! Move assignment operator
     ProjectionBase& operator=(ProjectionBase &&other) noexcept = default;
 
     //! initialises the fft engine (plan the transform)
     void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate);
 
     //! apply the projection operator to a field
     void apply_projection(Field_t & field) const;
 
+    //!
+    const Ccoord & get_resolutions() const {
+      return this->fft_engine.get_resolutions();}
+    const Rcoord & get_lengths() const {
+      return this->fft_engine.get_lengths();}
+
   protected:
     FFT_Engine & fft_engine;
     LFieldCollection_t & projection_container{};
 
   private:
   };
 
 }  // muSpectre
 
 
 
 #endif /* PROJECTION_BASE_H */
diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh
index 1da34a2..7f55cc3 100644
--- a/src/materials/material_base.hh
+++ b/src/materials/material_base.hh
@@ -1,139 +1,143 @@
 /**
  * file   material_base.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   25 Oct 2017
  *
  * @brief  Base class for materials (constitutive models)
  *
  * @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 <string>
 
 #include "common/common.hh"
 #include "common/field_map_matrixlike.hh"
 #include "common/field_collection.hh"
 
 
 #ifndef MATERIAL_BASE_H
 #define MATERIAL_BASE_H
 
 namespace muSpectre {
 
   //! DimS spatial dimension (dimension of problem
   //! DimM material_dimension (dimension of constitutive law)
   template<Dim_t DimS, Dim_t DimM>
   class MaterialBase
   {
   public:
     //! typedefs for data handled by this interface
     //! global field collection for system-wide fields, like stress, strain, etc
     using GFieldCollection_t = FieldCollection<DimS, DimM, true>;
     //! field collection for internal variables, such as eigen-strains,
     //! plastic strains, damage variables, etc, but also for managing which
     //! pixels the material is responsible for
     using MFieldCollection_t = FieldCollection<DimS, DimM, false>;
+    using iterator = typename MFieldCollection_t::iterator;
     using Field_t = internal::FieldBase<GFieldCollection_t>;
     using StressField_t = TensorField<GFieldCollection_t, Real, secondOrder, DimM>;
     using StrainField_t = StressField_t;
     using TangentField_t = TensorField<GFieldCollection_t, Real, fourthOrder, DimM>;
     using Ccoord = Ccoord_t<DimS>;
     //! Default constructor
     MaterialBase() = delete;
 
     //! Construct by name
     MaterialBase(std::string name);
 
     //! Copy constructor
     MaterialBase(const MaterialBase &other) = delete;
 
     //! Move constructor
     MaterialBase(MaterialBase &&other) noexcept = delete;
 
     //! Destructor
     virtual ~MaterialBase() noexcept = default;
 
     //! Copy assignment operator
     MaterialBase& operator=(const MaterialBase &other) = delete;
 
     //! Move assignment operator
     MaterialBase& operator=(MaterialBase &&other) noexcept = delete;
 
 
     /**
      *  take responsibility for a pixel identified by its cell coordinates
      *  WARNING: this won't work for materials with additional info per pixel
      *  (as, e.g. for eigenstrain), we need to pass more parameters. Materials
      *  of this tye need to overload add_pixel
      */
     void add_pixel(const Ccoord & ccooord);
 
     //! allocate memory, etc
     virtual void initialise(bool stiffness = false) = 0;
 
     //! return the materil's name
     const std::string & get_name() const;
 
     //! for static inheritance stuff
     constexpr static Dim_t sdim() {return DimS;}
     constexpr static Dim_t mdim() {return DimM;}
     //! computes stress
     virtual void compute_stresses(const StrainField_t & F,
                                   StressField_t & P,
                                   Formulation form) = 0;
     /**
      * Convenience function to compute stresses, mostly for debugging and
      * testing. Has runtime-cost associated with compatibility-checking and
      * conversion of the Field_t arguments that can be avoided by using the
      * version with strongly typed field references
      */
     void compute_stresses(const Field_t & F,
                           Field_t & P,
                           Formulation form);
     //! computes stress and tangent moduli
     virtual void compute_stresses_tangent(const StrainField_t & F,
                                           StressField_t & P,
                                           TangentField_t & K,
                                           Formulation form) = 0;
     /**
      * Convenience function to compute stresses and tangent moduli, mostly for
      * debugging and testing. Has runtime-cost associated with
      * compatibility-checking and conversion of the Field_t arguments that can
      * be avoided by using the version with strongly typed field references
      */
     void compute_stresses_tangent(const Field_t & F,
                                   Field_t & P,
                                   Field_t & K,
                                   Formulation form);
 
+    inline iterator begin() {return this->internal_fields.begin();}
+    inline iterator end()  {return this->internal_fields.end();}
+
   protected:
 
     //! members
     const std::string name;
     MFieldCollection_t internal_fields{};
 
 
   private:
   };
 }  // muSpectre
 
 #endif /* MATERIAL_BASE_H */
diff --git a/src/system/CMakeLists.txt b/src/system/CMakeLists.txt
index 3ff0f3c..a3c3189 100644
--- a/src/system/CMakeLists.txt
+++ b/src/system/CMakeLists.txt
@@ -1,9 +1,7 @@
 set (systems_SRC
   system_base.cc
-  #fftw_engine_c2c.cc
-  #fftw_engine_r2c.cc
   )
 find_package(FFTW REQUIRED)
 
 add_library(systems ${systems_SRC})
-target_link_libraries(systems Eigen3::Eigen common ${FFTW_LIBRARIES})
+target_link_libraries(systems Eigen3::Eigen common)
diff --git a/src/system/system_base.cc b/src/system/system_base.cc
index 8286acf..4051497 100644
--- a/src/system/system_base.cc
+++ b/src/system/system_base.cc
@@ -1,37 +1,127 @@
 /**
  * file   system_base.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   01 Nov 2017
  *
  * @brief  Implementation for system base class
  *
  * @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 <sstream>
+#include <algorithm>
+
 #include "system/system_base.hh"
+#include "common/ccoord_operations.hh"
 
 
 namespace muSpectre {
 
-  
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  SystemBase<DimS, DimM>::SystemBase(Projection_ptr projection)
+    :resolutions{projection->get_resolutions()},
+     lengths{projection->get_lengths()},
+     F{make_field<StrainField_t>("Gradient", this->fields)},
+     P{make_field<StressField_t>("Piola-Kirchhoff-1", this->fields)},
+     K_ptr{nullptr}, projection{std::move(projection)}
+  {}
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  void SystemBase<DimS, DimM>::add_material(Material_ptr mat) {
+    this->materials.push_back(std::move(mat));
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  typename SystemBase<DimS, DimM>::FullResponse_t
+  SystemBase<DimS, DimM>::evaluate_stress_tangent(StrainField_t & grad) {
+    if (this->is_initialised == false) {
+      this->initialise();
+    }
+    //! High level compatibility checks
+    if (grad.size() != this->F.size()) {
+      throw std::runtime_error("Size mismatch");
+    }
+    if (this->K_ptr == nullptr) {
+      K_ptr = &make_field<TangentField_t>("Tangent Stiffness", this->fields);
+    }
+
+    for (auto & mat: this->materials) {
+      mat->compute_stresses_tangent(grad, this->P, *this->K_ptr, form);
+    }
+    return std::tie(this->P, *this->K_ptr);
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  void SystemBase<DimS, DimM>::initialise() {
+    this->check_material_coverage();
+    this->fields.initialise(this->resolutions);
+    this->is_initialised = true;
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  void SystemBase<DimS, DimM>::check_material_coverage() {
+    auto nb_pixels = CcoordOps::get_size(this->resolutions);
+    std::vector<MaterialBase<DimS, DimM>*> assignments(nb_pixels, nullptr);
+    for (auto & mat: this->materials) {
+      for (auto & pixel: *mat) {
+        auto index = CcoordOps::get_index(this->resolutions, pixel);
+        auto& assignment{assignments.at(index)};
+        if (assignment != nullptr) {
+          std::stringstream err{};
+          err << "Pixel " << pixel << "is already assigned to material '"
+              << assignment->get_name()
+              << "' and cannot be reassigned to material '" << mat->get_name();
+          throw std::runtime_error(err.str());
+        } else {
+          assignments[index] = mat.get();
+        }
+      }
+    }
+
+    // find and identify unassigned pixels
+    std::vector<Ccoord> unassigned_pixels;
+    for (size_t i = 0; i < assignments.size(); ++i) {
+      if (assignments[i] == nullptr) {
+        unassigned_pixels.push_back(CcoordOps::get_ccoord(this->resolutions, i));
+      }
+    }
+
+    if (unassigned_pixels.size() != 0) {
+      std::stringstream err {};
+      err << "The following pixels have were not assigned a material: ";
+      for (auto & pixel: unassigned_pixels) {
+        err << pixel << ", ";
+      }
+      err << "and that cannot be handled";
+      throw std::runtime_error(err.str());
+    }
+  }
+
+  template class SystemBase<twoD, twoD>;
+  template class SystemBase<threeD, threeD>;
 
 }  // muSpectre
diff --git a/src/system/system_base.hh b/src/system/system_base.hh
index d6afc2c..d35adc9 100644
--- a/src/system/system_base.hh
+++ b/src/system/system_base.hh
@@ -1,46 +1,113 @@
 /**
  * file   system_base.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   01 Nov 2017
  *
- * @brief  Base class representing a a unit cell system
+ * @brief Base class representing a unit cell system with single
+ *        projection operator
  *
  * @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 <common/common.hh>
+#include <vector>
+#include <memory>
+#include <tuple>
+
+#include "common/common.hh"
+#include "common/field_collection.hh"
+#include "materials/material_base.hh"
+#include "fft/projection_base.hh"
+
 #ifndef SYSTEM_BASE_H
 #define SYSTEM_BASE_H
 
 namespace muSpectre {
 
   //! DimS spatial dimension (dimension of problem
   //! DimM material_dimension (dimension of constitutive law)
   template<Dim_t DimS, Dim_t DimM>
   class SystemBase
   {
+  public:
+    constexpr static Formulation form{Formulation::finite_strain};
+    using Ccoord = Ccoord_t<DimS>;
+    using Rcoord = Rcoord_t<DimS>;
+    using FieldCollection_t = FieldCollection<DimS, DimM>;
+    using Material_ptr = std::unique_ptr<MaterialBase<DimS, DimM>>;
+    using Projection_ptr = std::unique_ptr<ProjectionBase<DimS, DimM>>;
+    using StrainField_t = TensorField<FieldCollection_t, Real, secondOrder, DimM>;
+    using StressField_t = TensorField<FieldCollection_t, Real, secondOrder, DimM>;
+    using TangentField_t = TensorField<FieldCollection_t, Real, secondOrder, DimM>;
+    using FullResponse_t = std::tuple<const StressField_t&, const TangentField_t&>;
+
+    //! Default constructor
+    SystemBase() = delete;
+
+    //! constructor using sizes and resolution
+    SystemBase(Projection_ptr projection);
+
+    //! Copy constructor
+    SystemBase(const SystemBase &other) = delete;
+
+    //! Move constructor
+    SystemBase(SystemBase &&other) noexcept = default;
+
+    //! Destructor
+    virtual ~SystemBase() noexcept = default;
+
+    //! Copy assignment operator
+    SystemBase& operator=(const SystemBase &other) = delete;
+
+    //! Move assignment operator
+    SystemBase& operator=(SystemBase &&other) noexcept = default;
+
+    /**
+     * Materials can only be moved. This is to assure exclusive
+     * ownership of any material by this system
+     */
+    void add_material(Material_ptr mat);
+
+    FullResponse_t evaluate_stress_tangent(StrainField_t & F);
+
+    void initialise();
+  protected:
+    //! make sure that every pixel is assigned to one and only one material
+    void check_material_coverage();
+
+    const Ccoord & resolutions;
+    const Rcoord & lengths;
+    FieldCollection_t fields{};
+    StrainField_t & F;
+    StressField_t & P;
+    //! Tangent field migth not even be required; so this is a
+    //! pointer instead of a ref
+    TangentField_t * K_ptr;
+    std::vector<Material_ptr> materials{};
+    Projection_ptr projection;
+    bool is_initialised{false};
+  private:
   };
 
 
 }  // muSpectre
 
 #endif /* SYSTEM_BASE_H */
diff --git a/tests/test_system_base.cc b/tests/test_system_base.cc
new file mode 100644
index 0000000..8870c73
--- /dev/null
+++ b/tests/test_system_base.cc
@@ -0,0 +1,80 @@
+/**
+ * file   test_system_base.cc
+ *
+ * @author Till Junge <till.junge@epfl.ch>
+ *
+ * @date   14 Dec 2017
+ *
+ * @brief  Tests for the basic system class
+ *
+ * @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 <boost/mpl/list.hpp>
+
+#include "tests.hh"
+#include "common/common.hh"
+#include "system/system_base.hh"
+#include "fft/projection_finite_strain_fast.hh"
+#include "fft/fftw_engine.hh"
+
+
+namespace muSpectre {
+
+  BOOST_AUTO_TEST_SUITE(system_base);
+  template <Dim_t DimS>
+  struct Sizes {
+  };
+  template<>
+  struct Sizes<twoD> {
+    constexpr static Ccoord_t<twoD> get_resolution() {
+      return Ccoord_t<twoD>{3, 5};}
+    constexpr static Rcoord_t<twoD> get_lengths() {
+      return Rcoord_t<twoD>{3.4, 5.8};}
+  };
+  template<>
+  struct Sizes<threeD> {
+    constexpr static Ccoord_t<threeD> get_resolution() {
+      return Ccoord_t<threeD>{3, 5, 7};}
+    constexpr static Rcoord_t<threeD> get_lengths() {
+      return Rcoord_t<threeD>{3.4, 5.8, 6.7};}
+  };
+
+  template <Dim_t DimS, Dim_t DimM>
+  struct SystemBaseFixture: SystemBase<DimS, DimM> {
+    SystemBaseFixture()
+      :fft_engine{Sizes<DimS>::get_resolution(), Sizes<DimS>::get_lengths()},
+       proj{fft_engine},
+       SystemBase<DimS, DimM>{std::move(proj)}{}
+
+    FFTW_Engine<DimS, DimM> fft_engine;
+    ProjectionFiniteStrainFast<DimS, DimM> proj;
+  };
+
+  using fixlist = boost::mpl::list<SystemBaseFixture<twoD, twoD>,
+                                   SystemBaseFixture<threeD, threeD>>;
+
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) {
+  }
+
+  BOOST_AUTO_TEST_SUITE_END();
+
+}  // muSpectre