diff --git a/src/common/statefield.hh b/src/common/statefield.hh
index 19b11ff..1337cce 100644
--- a/src/common/statefield.hh
+++ b/src/common/statefield.hh
@@ -1,546 +1,547 @@
 /**
  * file   statefield.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   28 Feb 2018
  *
  * @brief A state field is an abstraction of a field that can hold
  * current, as well as a chosen number of previous values. This is
  * useful for instance for internal state variables in plastic laws,
  * where a current, new, or trial state is computed based on its
  * previous state, and at convergence, this new state gets cycled into
  * the old, the old into the old-1 etc. The state field abstraction
  * helps doing this safely (i.e. only const references to the old
  * states are available, while the current state can be assigned
  * to/modified), and efficiently (i.e., no need to copy values from
  * new to old, we just cycle the labels). This file implements the
  * state field as well as state maps using the Field, FieldCollection
  * and FieldMap abstractions of µSpectre
  *
  * Copyright © 2018 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 STATEFIELD_H
 #define STATEFIELD_H
 
 #include "common/field.hh"
 #include "common/utilities.hh"
 
 #include <array>
 #include <string>
 #include <sstream>
 
 namespace muSpectre {
 
   /**
    * Base class for state fields, useful for storing polymorphic references
    */
   template <size_t nb_memory>
   class StateFieldBase {
   public:
     //! the current (historically accurate) ordering of the fields
     using index_t = std::array<size_t, nb_memory+1>;
     //! get the current ordering of the fields
     inline const index_t & get_indices() const {return this->indices;}
   protected:
     //! constructor
     StateFieldBase(index_t indices):indices{indices}{};
     index_t indices; ///< these are cycled through
   };
 
   //! early declaration
   template <class FieldMap, size_t nb_memory>
   class StateFieldMap;
 
   namespace internal {
 
     template <class Field, size_t size, size_t... I>
     inline decltype(auto)
     build_fields_helper(std::string prefix,
                         typename Field::Base::collection_t & collection,
                         std::index_sequence<I...>) {
       auto get_field{[&prefix, &collection](size_t i) -> Field&
           {
             std::stringstream name_stream{};
             name_stream << prefix << ", sub_field index " << i;
             return make_field<Field>(name_stream.str(), collection);
           }};
       return std::tie(get_field(I)...);
     }
 
     /* ---------------------------------------------------------------------- */
     template <size_t size, size_t... I>
     decltype(auto) build_indices(std::index_sequence<I...>) {
       return std::array<size_t, size>{I...};
     }
 
   }  // internal
 
   /**
    * A statefield is an abstraction around a Field that can hold a
    * current and `nb_memory` previous values. There are useful for
    * history variables, for instance.
    */
   template <class Field_t, size_t nb_memory=1>
   class StateField: public StateFieldBase<nb_memory>
   {
   public:
     //! Base class for all state fields of same memory 
     using Base = StateFieldBase<nb_memory>;
     //! the underlying field's collection type
     using FieldCollection = typename Field_t::Base::collection_t;
     /**
      * storage of field refs (can't be a `std::array`, because arrays
      * of refs are explicitely forbidden
      */
     using Fields_t = tuple_array<Field_t&, nb_memory+1>;
 
     //! Default constructor
     StateField() = delete;
 
     /**
      * Constructor.  @param unique_prefix is used to create the names
      * of the fields that this abstraction creates in the background
      * @param collection is the field collection in which the
      * subfields will be stored
      */
     inline StateField(std::string unique_prefix, FieldCollection & collection)
       : Base{internal::build_indices<nb_memory+1>(std::make_index_sequence<nb_memory+1>{})},
         prefix{unique_prefix}, fields{internal::build_fields_helper<Field_t, nb_memory+1>
         (unique_prefix, collection, std::make_index_sequence<nb_memory+1>{})}
     {}
 
     //! Copy constructor
     StateField(const StateField &other) = delete;
 
-    //! Move constructor
+    //! Move coonstructor
     StateField(StateField &&other) = delete;
 
     //! Destructor
     virtual ~StateField() = default;
 
     //! Copy assignment operator
     StateField& operator=(const StateField &other) = delete;
 
     //! Move assignment operator
     StateField& operator=(StateField &&other) = delete;
 
     //! get (modifiable) current field
-    inline StateField& current() {
+    inline Field_t& current() {
+      // TODO: doesn't work
       return this->fields[this->indices[0]];
     }
 
     //! get (constant) previous field
     template <size_t nb_steps_ago=1>
-    inline const StateField& old() {
+    inline const Field_t& old() {
       static_assert(nb_steps_ago <= nb_memory,
                     "you can't go that far inte the past");
       static_assert(nb_steps_ago > 0,
                     "Did you mean to call current()?");
       return this->fields[this->indices[nb_steps_ago]];
     }
 
     // //! factory function
     // template<class StateFieldType, class CollectionType, typename... Args>
     // friend StateFieldType& make_statefield(std::string unique_prefix,
     //                                        CollectionType & collection,
     //                                        Args&&... args);
 
     //! returns a `StateField` reference if `other is a compatible state field
     inline static StateField& check_ref(Base& other) {
       // the following triggers and exception if the fields are incompatible
       Field_t::check_ref(other.fields[0]);
       return static_cast<StateField&> (other);
     }
 
     //! returns a const `StateField` reference if `other` is a compatible state field
     inline static const StateField& check_ref(const Base& other) {
       // the following triggers and exception if the fields are incompatible
       Field_t::check_ref(other.fields[0]);
       return static_cast<const StateField&> (other);
     }
 
     //! get naming prefix
     const std::string & get_prefix() const {return this->prefix;}
 
     //! get a ref to the `StateField` 's field collection
     const FieldCollection & get_collection() const {
       return std::get<0>(this->fields).get_collection();}
 
     //! get a ref to the `StateField` 's fields
      Fields_t & get_fields() {
        return this->fields;
      }
 
     /**
      * Pure convenience functions to get a MatrixFieldMap of
      * appropriate dimensions mapped to this field. You can also
      * create other types of maps, as long as they have the right
      * fundamental type (T), the correct size (nbComponents), and
      * memory (nb_memory).
      */
     inline decltype(auto) get_map() {
       using FieldMap = decltype(std::get<0>(this->fields).get_map());
       return StateFieldMap<FieldMap, nb_memory>(*this);
     }
 
     /**
      * Pure convenience functions to get a MatrixFieldMap of
      * appropriate dimensions mapped to this field. You can also
      * create other types of maps, as long as they have the right
      * fundamental type (T), the correct size (nbComponents), and
      * memory (nb_memory).
      */
     inline decltype(auto) get_const_map() {
       using FieldMap = decltype(std::get<0>(this->fields).get_const_map());
       return StateFieldMap<FieldMap, nb_memory>(*this);
     }
 
     /**
      * cycle the fields (current becomes old, old becomes older,
      * oldest becomes current)
      */
     inline void cycle() {
       for (auto & val: this->indices) {
         val = (val+1)%(nb_memory+1);
       }
     }
 
   protected:
     /**
      * the unique prefix is used as the first part of the unique name
      * of the subfields belonging to this state field
      */
     std::string prefix;
     Fields_t fields; //!< container for the states
   private:
   };
 
 
   namespace internal {
 
     template <class FieldMap, size_t size, class Fields, size_t... I>
     inline decltype(auto) build_maps_helper(Fields & fields,
                                             std::index_sequence<I...>) {
       return std::array<FieldMap, size>{FieldMap(std::get<I>(fields))...};
     }
 
   }  // internal
 
   /**
    * extends the StateField <-> Field equivalence to StateFieldMap <-> FieldMap
    */
   template <class FieldMap, size_t nb_memory=1>
   class StateFieldMap
   {
   public:
     /**
      * iterates over all pixels in the `muSpectre::FieldCollection` and
      * dereferences to a proxy giving access to the appropriate iterates
      * of the underlying `FieldMap` type.
      */
     class iterator;
 
     //! stl conformance
     using reference = typename iterator::reference;
     //! stl conformance
     using value_type = typename iterator::value_type;
     //! stl conformance
     using size_type = typename iterator::size_type;
 
     //! field collection type where this state field can be stored
     using FieldCollection= typename FieldMap::Field::collection_t;
     //! base class (for polymorphic references)
     using StateFieldBase_t = StateFieldBase<nb_memory>;
     //! for traits access
     using FieldMap_t = FieldMap;
     //! for traits access
     using ConstFieldMap_t = typename FieldMap::ConstMap;
 
     //! Default constructor
     StateFieldMap() = delete;
 
     //! constructor using a StateField
     template <class StateField>
     StateFieldMap(StateField & statefield)
       :collection{statefield.get_collection()},
        statefield{statefield},
        maps{internal::build_maps_helper
            <FieldMap, nb_memory+1>(statefield.get_fields(),
                                    std::make_index_sequence<nb_memory+1>{})},
        const_maps{internal::build_maps_helper
            <ConstFieldMap_t, nb_memory+1>(statefield.get_fields(),
                                           std::make_index_sequence<nb_memory+1>{})}
     {
       static_assert(std::is_base_of<StateFieldBase_t, StateField>::value,
                     "Not the right type of StateField ref");
     }
 
     //! Copy constructor
     StateFieldMap(const StateFieldMap &other) = delete;
 
     //! Move constructor
     StateFieldMap(StateFieldMap &&other) = default;
 
     //! Destructor
     virtual ~StateFieldMap() = default;
 
     //! Copy assignment operator
     StateFieldMap& operator=(const StateFieldMap &other) = delete;
 
     //! Move assignment operator
     StateFieldMap& operator=(StateFieldMap &&other) = delete;
 
     //! access the wrapper to a given pixel directly
     value_type operator[](size_type index) {
       return *iterator(*this, index);
     }
 
     /**
      * return a ref to the current field map. useful for instance for
      * initialisations of `StateField` instances
      */
     FieldMap& current() {
       return this->maps[this->statefield.get_indices()[0]];
     }
 
     //! stl conformance
     iterator begin() {
       return iterator(*this, 0);}
     //! stl conformance
     iterator end() {
       return iterator(*this, this->collection.size());}
 
   protected:
     const FieldCollection & collection; //!< collection holding the field
     StateFieldBase_t & statefield; //!< ref to the field itself
     std::array<FieldMap, nb_memory+1> maps;//!< refs to the addressable maps;
     //! const refs to the addressable maps;
     std::array<ConstFieldMap_t, nb_memory+1> const_maps;
   private:
   };
 
   /**
    * Iterator class used by the `StateFieldMap`
    */
   template <class FieldMap, size_t nb_memory>
   class StateFieldMap<FieldMap, nb_memory>::iterator
   {
   public:
     class StateWrapper;
 
     using Ccoord = typename FieldMap::Ccoord; //!< cell coordinates type
     using value_type = StateWrapper; //!< stl conformance
     using const_value_type = value_type; //!< stl conformance
     using pointer_type = value_type*; //!< stl conformance
     using difference_type = std::ptrdiff_t; //!< stl conformance
     using size_type = size_t; //!< stl conformance
     using iterator_category = std::random_access_iterator_tag; //!< stl conformance
     using reference = StateWrapper; //!< stl conformance
 
     //! Default constructor
     iterator() = delete;
 
     //! constructor
     iterator(StateFieldMap& map, size_t index = 0)
       :index{index}, map{map}
     {};
 
     //! Copy constructor
     iterator(const iterator &other) = default;
 
     //! Move constructor
     iterator(iterator &&other) = default;
 
     //! Destructor
     virtual ~iterator() = default;
 
     //! Copy assignment operator
     iterator& operator=(const iterator &other) = default;
 
     //! Move assignment operator
     iterator& operator=(iterator &&other) = default;
 
     //! pre-increment
     inline iterator & operator++() {
       this->index++; return *this;}
     //! post-increment
     inline iterator operator++(int) {
       iterator curr{*this}; this->index++; return curr;}
     //! dereference
     inline value_type operator*() {
       return value_type(*this);}
     //! pre-decrement
     inline iterator & operator--() {
       this->index--; return *this;}
     //! post-decrement
     inline iterator operator--(int) {
       iterator curr{*this}; this->index--; return curr;}
     //! access subscripting
     inline value_type operator[](difference_type diff) {
       return value_type{iterator{this->map, this->index+diff}};}
     //! equality
     inline bool operator==(const iterator & other) const {
       return this->index == other.index;
     }
     //! inequality
     inline bool operator!=(const iterator & other) const {
       return this->index != other.index;}
     //! div. comparisons
     inline bool operator<(const iterator & other) const {
       return this->index < other.index;
     }
     //! div. comparisons
     inline bool operator<=(const iterator & other) const {
       return this->index <= other.index;
     }
     //! div. comparisons
     inline bool operator>(const iterator & other) const {
       return this->index > other.index;
     }
     //! div. comparisons
     inline bool operator>=(const iterator & other) const {
       return this->index >= other.index;
     }
     //! additions, subtractions and corresponding assignments
     inline iterator operator+(difference_type diff) const {
       return iterator{this->map, this-index + diff};
     }
     //! additions, subtractions and corresponding assignments
     inline iterator operator-(difference_type diff) const {
       return iterator{this->map, this-index - diff};}
     //! additions, subtractions and corresponding assignments
     inline iterator& operator+=(difference_type diff) {
       this->index += diff; return *this;}
     //! additions, subtractions and corresponding assignments
     inline iterator& operator-=(difference_type diff) {
       this->index -= diff; return *this;
     }
 
     //! get pixel coordinates
     inline Ccoord get_ccoord() const {
       return this->map.collection.get_ccoord(this->index);
     }
 
     //! access the index
     inline const size_t & get_index() const {return this->index;}
 
   protected:
     size_t index; //!< current pixel this iterator refers to
     StateFieldMap& map; //!< map over with `this` iterates
   private:
   };
 
 
 
   namespace internal {
 
     //! FieldMap is an `Eigen::Map` or `Eigen::TensorMap` here
     template <class FieldMap, size_t size, size_t... I,
               class iterator, class maps_t, class indices_t>
     inline decltype(auto)
     build_old_vals_helper(iterator& it, maps_t & maps, indices_t & indices,
                           std::index_sequence<I...>) {
       return tuple_array<FieldMap, size>(std::forward_as_tuple(maps[indices[I+1]][it.get_index()]...));
     }
 
     template <class FieldMap, size_t size, class iterator, class maps_t, class indices_t>
     inline decltype(auto)
     build_old_vals(iterator& it, maps_t & maps, indices_t & indices) {
       return tuple_array<FieldMap, size>{build_old_vals_helper<FieldMap, size>
           (it, maps, indices, std::make_index_sequence<size>{})};
     }
 
   }  // internal
 
   /**
    * Light-weight resource-handle representing the current and old
    * values of a field at a given pixel identified by an iterator
    * pointing to it
    */
   template <class FieldMap, size_t nb_memory>
   class StateFieldMap<FieldMap, nb_memory>::iterator::StateWrapper
   {
   public:
     //! short-hand
     using iterator = typename StateFieldMap::iterator;
     //! short-hand
     using Ccoord = typename iterator::Ccoord;
     //! short-hand
     using Map = typename FieldMap::reference;
     //! short-hand
     using ConstMap = typename FieldMap::const_reference;
 
     //! Default constructor
     StateWrapper() = delete;
 
     //! Copy constructor
     StateWrapper(const StateWrapper &other) = default;
 
     //! Move constructor
     StateWrapper(StateWrapper &&other) = default;
 
     //! construct with `StateFieldMap::iterator`
     StateWrapper(iterator & it)
       :it{it},
        current_val{it.map.maps[it.map.statefield.get_indices()[0]][it.index]},
        old_vals(internal::build_old_vals
            <ConstMap, nb_memory>(it, it.map.const_maps,
                                  it.map.statefield.get_indices()))
     {    }
 
     //! Destructor
     virtual ~StateWrapper() = default;
 
     //! Copy assignment operator
     StateWrapper& operator=(const StateWrapper &other) = default;
 
     //! Move assignment operator
     StateWrapper& operator=(StateWrapper &&other) = default;
 
     //! returns reference to the currectly mapped value
     inline Map& current() {
       return this->current_val;
     }
 
     //! recurnts reference the the value that was current `nb_steps_ago` ago
     template <size_t nb_steps_ago = 1>
     inline const ConstMap & old() const{
       static_assert (nb_steps_ago <= nb_memory,
                      "You have not stored that time step");
       static_assert (nb_steps_ago > 0,
                      "Did you mean to access the current value? If so, use "
                      "current()");
       return std::get<nb_memory-nb_steps_ago>(this->old_vals);
     }
 
     //! read the coordinates of the current pixel
     inline Ccoord get_ccoord() const {
       return this->it.get_ccoord();
     }
 
   protected:
     iterator& it; //!< ref to the iterator that dereferences to `this`
     Map current_val; //!< current value
     tuple_array<ConstMap, nb_memory> old_vals; //!< all stored old values
   private:
   };
 }  // muSpectre
 
 #endif /* STATEFIELD_H */
diff --git a/src/common/tensor_algebra.hh b/src/common/tensor_algebra.hh
index 4439c39..0de523d 100644
--- a/src/common/tensor_algebra.hh
+++ b/src/common/tensor_algebra.hh
@@ -1,342 +1,351 @@
 /**
  * @file   tensor_algebra.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   05 Nov 2017
  *
  * @brief  collection of compile-time quantities and algrebraic functions for
  *         tensor operations
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #ifndef TENSOR_ALGEBRA_H
 #define TENSOR_ALGEBRA_H
 
 #include "common/T4_map_proxy.hh"
 #include "common/common.hh"
 #include "common/eigen_tools.hh"
 
 #include <Eigen/Dense>
 #include <unsupported/Eigen/CXX11/Tensor>
 
 #include <type_traits>
 
 namespace muSpectre {
 
   namespace Tensors {
 
     //! second-order tensor representation
     template<Dim_t dim>
     using Tens2_t = Eigen::TensorFixedSize<Real, Eigen::Sizes<dim, dim>>;
     //! fourth-order tensor representation
     template<Dim_t dim>
     using Tens4_t = Eigen::TensorFixedSize<Real, Eigen::Sizes<dim, dim, dim, dim>>;
 
     //----------------------------------------------------------------------------//
     //! compile-time second-order identity
     template<Dim_t dim>
     constexpr inline Tens2_t<dim> I2() {
       Tens2_t<dim> T;
       using Mat_t = Eigen::Matrix<Real, dim, dim>;
       Eigen::Map<Mat_t>(&T(0, 0)) = Mat_t::Identity();
       return T;
     }
 
     /* ---------------------------------------------------------------------- */
     //! Check whether a given expression represents a Tensor specified order
     template<class T, Dim_t order>
     struct is_tensor {
       //! evaluated test
       constexpr static bool value = (std::is_convertible
                       <T, Eigen::Tensor<Real, order>>::value ||
                     std::is_convertible
                       <T, Eigen::Tensor<Int, order>>::value ||
                     std::is_convertible
                       <T, Eigen::Tensor<Complex, order>>::value);
     };
     /* ---------------------------------------------------------------------- */
     /** compile-time outer tensor product as defined by Curnier
      *  R_ijkl = A_ij.B_klxx
      *    0123     01   23
      */
     template<Dim_t dim, typename T1, typename T2>
     constexpr inline decltype(auto) outer(T1 && A, T2 && B) {
       // Just make sure that the right type of parameters have been given
       constexpr Dim_t order{2};
       static_assert(is_tensor<T1, order>::value,
                     "T1 needs to be convertible to a second order Tensor");
       static_assert(is_tensor<T2, order>::value,
                     "T2 needs to be convertible to a second order Tensor");
 
       // actual function
       std::array<Eigen::IndexPair<Dim_t>, 0> dims{};
       return A.contract(B, dims);
     }
 
     /* ---------------------------------------------------------------------- */
     /** compile-time underlined outer tensor product as defined by Curnier
      *  R_ijkl = A_ik.B_jlxx
      *    0123     02   13
      *    0213     01   23 <- this defines the shuffle order
      */
     template<Dim_t dim, typename T1, typename T2>
     constexpr inline decltype(auto) outer_under(T1 && A, T2 && B) {
       constexpr size_t order{4};
       return outer<dim>(A, B).shuffle(std::array<Dim_t, order>{{0, 2, 1, 3}});
     }
 
     /* ---------------------------------------------------------------------- */
     /** compile-time overlined outer tensor product as defined by Curnier
      *  R_ijkl = A_il.B_jkxx
      *    0123     03   12
      *    0231     01   23 <- this defines the shuffle order
      */
     template<Dim_t dim, typename T1, typename T2>
     constexpr inline decltype(auto) outer_over(T1 && A, T2 && B) {
       constexpr size_t order{4};
       return outer<dim>(A, B).shuffle(std::array<Dim_t, order>{{0, 2, 3, 1}});
     }
 
     //! compile-time fourth-order symmetrising identity
     template<Dim_t dim>
     constexpr inline Tens4_t<dim> I4S() {
       auto I = I2<dim>();
       return 0.5*(outer_under<dim>(I, I) + outer_over<dim>(I, I));
     }
 
   }  // Tensors
 
 
   namespace Matrices {
 
     //! second-order tensor representation
     template<Dim_t dim>
     using Tens2_t = Eigen::Matrix<Real, dim, dim>;
     //! fourth-order tensor representation
     template<Dim_t dim>
     using Tens4_t = T4Mat<Real, dim>;
 
     //----------------------------------------------------------------------------//
     //! compile-time second-order identity
     template<Dim_t dim>
     constexpr inline Tens2_t<dim> I2() {
       return Tens2_t<dim>::Identity();
     }
 
     /* ---------------------------------------------------------------------- */
     /** compile-time outer tensor product as defined by Curnier
      *  R_ijkl = A_ij.B_klxx
      *    0123     01   23
      */
     template<typename T1, typename T2>
     constexpr inline decltype(auto) outer(T1 && A, T2 && B) {
       // Just make sure that the right type of parameters have been given
       constexpr Dim_t dim{EigenCheck::tensor_dim<T1>::value};
       static_assert((dim == EigenCheck::tensor_dim<T2>::value),
                      "A and B do not have the same dimension");
 
       Tens4_t<dim> product;
 
       for (Dim_t i = 0; i < dim; ++i) {
         for (Dim_t j = 0; j < dim; ++j) {
           for (Dim_t k = 0; k < dim; ++k) {
             for (Dim_t l = 0; l < dim; ++l) {
               get(product, i, j, k, l) = A(i, j) * B(k, l);
             }
           }
         }
       }
       return product;
     }
 
     /* ---------------------------------------------------------------------- */
     /** compile-time underlined outer tensor product as defined by Curnier
      *  R_ijkl = A_ik.B_jlxx
      *    0123     02   13
      *    0213     01   23 <- this defines the shuffle order
      */
     template<typename T1, typename T2>
     constexpr inline decltype(auto) outer_under(T1 && A, T2 && B) {
       // Just make sure that the right type of parameters have been given
       constexpr Dim_t dim{EigenCheck::tensor_dim<T1>::value};
       static_assert((dim == EigenCheck::tensor_dim<T2>::value),
                      "A and B do not have the same dimension");
 
       Tens4_t<dim> product;
 
       for (Dim_t i = 0; i < dim; ++i) {
         for (Dim_t j = 0; j < dim; ++j) {
           for (Dim_t k = 0; k < dim; ++k) {
             for (Dim_t l = 0; l < dim; ++l) {
               get(product, i, j, k, l) = A(i, k) * B(j, l);
             }
           }
         }
       }
       return product;
     }
 
     /* ---------------------------------------------------------------------- */
     /** compile-time overlined outer tensor product as defined by Curnier
      *  R_ijkl = A_il.B_jkxx
      *    0123     03   12
      *    0231     01   23 <- this defines the shuffle order
      */
     template<typename T1, typename T2>
     constexpr inline decltype(auto) outer_over(T1 && A, T2 && B) {
       // Just make sure that the right type of parameters have been given
       constexpr Dim_t dim{EigenCheck::tensor_dim<T1>::value};
       static_assert((dim == EigenCheck::tensor_dim<T2>::value),
                      "A and B do not have the same dimension");
 
       Tens4_t<dim> product;
 
       for (Dim_t i = 0; i < dim; ++i) {
         for (Dim_t j = 0; j < dim; ++j) {
           for (Dim_t k = 0; k < dim; ++k) {
             for (Dim_t l = 0; l < dim; ++l) {
               get(product, i, j, k, l) = A(i, l) * B(j, k);
             }
           }
         }
       }
       return product;
     }
 
 
     /**
      * Standart tensor multiplication
      */
     template<typename T4, typename T2>
     constexpr inline decltype(auto) tensmult(T4 && A, T2 && B) {
       constexpr Dim_t dim{EigenCheck::tensor_4_dim<T4>::value};
       static_assert((dim == EigenCheck::tensor_dim<T2>::value),
                     "Dimensionality check failed. Expects A to be a fourth-"
                     "order tensor of dimension dim and B to be a second-order "
                     "Tensor of same dimension");
       Tens2_t<dim> result;
       result.setZero();
 
       for (Dim_t i = 0; i < dim; ++i) {
         for (Dim_t j = 0; j < dim; ++j) {
           for (Dim_t k = 0; k < dim; ++k) {
             for (Dim_t l = 0; l < dim; ++l) {
               result(i,j) +=get(A, i, j, k, l) * B(k, l);
             }
           }
         }
       }
       return result;
     }
 
     //! compile-time fourth-order tracer
     template <Dim_t dim>
     constexpr inline Tens4_t<dim> Itrac() {
       auto I = I2<dim>();
       return outer(I,I);
     }
 
     //! compile-time fourth-order identity
     template <Dim_t dim>
     constexpr inline Tens4_t<dim> Iiden() {
       auto I = I2<dim>();
       return outer_under(I,I);
     }
 
     //! compile-time fourth-order transposer
     template <Dim_t dim>
     constexpr inline Tens4_t<dim> Itrns() {
       auto I = I2<dim>();
       return outer_over(I,I);
     }
 
     //! compile-time fourth-order symmetriser
     template<Dim_t dim>
     constexpr inline Tens4_t<dim> Isymm() {
       auto I = I2<dim>();
       return 0.5*(outer_under(I, I) + outer_over(I, I));
     }
 
     namespace internal {
 
       /* ---------------------------------------------------------------------- */
       template <Dim_t Dim, Dim_t Rank1, Dim_t Rank2>
       struct Dotter {
       };
 
       /* ---------------------------------------------------------------------- */
       template<Dim_t Dim>
       struct Dotter<Dim, secondOrder, fourthOrder> {
         template <class T1, class T2>
         static constexpr decltype(auto) dot(T1 && t1, T2 && t2) {
           using T4_t = T4Mat<typename std::remove_reference_t<T1>::Scalar, Dim>;
           T4_t ret_val{T4_t::Zero()};
           for (Int i = 0; i < Dim; ++i) {
             for (Int a = 0; a < Dim; ++a) {
               for (Int j = 0; j < Dim; ++j) {
                 for (Int k = 0; k < Dim; ++k) {
                   for (Int l = 0; l < Dim; ++l) {
                     for (Int m = 0; m < Dim; ++m) {
                       get(ret_val, i,j,k,l) +=
                         t1(i,a) * get(t2,a,j,k,l);
                     }
                   }
                 }
               }
             }
           }
           return ret_val;
         }
       };
 
       /* ---------------------------------------------------------------------- */
       template<Dim_t Dim>
       struct Dotter<Dim, fourthOrder, fourthOrder> {
         template <class T1, class T2>
         static constexpr decltype(auto) ddot(T1 && t1, T2 && t2) {
           return t1*t2;
         }
       };
 
+      /* ---------------------------------------------------------------------- */
+      template<Dim_t Dim>
+      struct Dotter<Dim, secondOrder, secondOrder> {
+        template <class T1, class T2>
+        static constexpr decltype(auto) ddot(T1 && t1, T2 && t2) {
+          return (t1*t2.transpose()).trace();
+        }
+      };
+
     }  // internal
 
     /* ---------------------------------------------------------------------- */
     template <Dim_t Dim, class T1, class T2>
     decltype(auto) dot(T1 && t1, T2 && t2) {
       constexpr Dim_t rank1{EigenCheck::tensor_rank<T1, Dim>::value};
       constexpr Dim_t rank2{EigenCheck::tensor_rank<T2, Dim>::value};
       return internal::Dotter<Dim, rank1, rank2>::dot(std::forward<T1>(t1),
                                                       std::forward<T2>(t2));
     }
 
     /* ---------------------------------------------------------------------- */
     template <Dim_t Dim, class T1, class T2>
     decltype(auto) ddot(T1 && t1, T2 && t2) {
       constexpr Dim_t rank1{EigenCheck::tensor_rank<T1, Dim>::value};
       constexpr Dim_t rank2{EigenCheck::tensor_rank<T2, Dim>::value};
       return internal::Dotter<Dim, rank1, rank2>::ddot(std::forward<T1>(t1),
                                                        std::forward<T2>(t2));
     }
 
   }  // Matrices
 }  // muSpectre
 
 
 #endif /* TENSOR_ALGEBRA_H */
diff --git a/src/materials/material_crystal_plasticity_finite.cc b/src/materials/material_crystal_plasticity_finite.cc
index 9177475..f1ab9ea 100644
--- a/src/materials/material_crystal_plasticity_finite.cc
+++ b/src/materials/material_crystal_plasticity_finite.cc
@@ -1,88 +1,129 @@
 /**
  * @file   material_crystal_plasticity_finite.cc
  *
  * @author Till Junge <till.junge@altermail.ch>
  * @author Francesco Maresca <francesco.maresca@epfl.ch>
  *
  * @date   23 Feb 2018
  *
  * @brief finite strain crystal plasticity implementation
  *
  * Copyright © 2018 Till Junge, Francesco Maresca
  *
  * µ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 "materials/material_crystal_plasticity_finite.hh"
+#include "common/iterators.hh"
 
 namespace muSpectre {
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM, Int NbSlip>
   MaterialCrystalPlasticityFinite<DimS, DimM, NbSlip>::
   MaterialCrystalPlasticityFinite(std::string name, Real bulk_m, Real shear_m, Real gamma_dot0, Real m_par, Real tau_y0, Real h0, Real s_infty, Real a_par, Real q_n, SlipVecs_ref Slip0, SlipVecs_ref Normal0, Real delta_t, Real tolerance, Int maxiter)
     : Parent{name},
       FpField("Plastic Deformation Gradient Fₚ(t)",this->internal_fields),
       GammaDotField("Plastic slip rates dγᵅ/dt",this->internal_fields),
       TauYField("Critical resolved shear stress τᵅy(t)",this->internal_fields),
       GammaField{
         make_field<MatrixField<LColl_t, Real, NbSlip, 1>>
           ("Accumulated slips γᵅ(t)",this->internal_fields)},
       EulerField{
         make_field<MatrixField<LColl_t, Real, NbEuler, 1>>
           ("Euler angles", this->internal_fields)},
       bulk_m{bulk_m}, shear_m{shear_m}, gamma_dot0{gamma_dot0}, m_par{m_par},
       tau_y0{tau_y0}, h0{h0}, s_infty{s_infty}, a_par{a_par}, q_n{q_n},
       delta_t{delta_t}, tolerance{tolerance}, maxiter{maxiter}, Slip0{Slip0},
       Normal0{Normal0},
       internal_variables{FpField.get_map(),
           GammaDotField.get_map(),
           TauYField.get_map(),
-          GammaField.get_map(),
           ArrayFieldMap<LColl_t, Real, NbEuler, 1, true>(EulerField)}
   {
     // Enforce n₀ and s₀ to be unit vectors!
     auto lambda{MatTB::convert_elastic_modulus<ElasticModulus::lambda,
                                                ElasticModulus::Bulk,
                                                ElasticModulus::Shear>(bulk_m,shear_m)};
     this->C_el=Hooke::compute_C_T4(lambda,shear_m);
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM, Int NbSlip>
   void MaterialCrystalPlasticityFinite<DimS, DimM, NbSlip>::
   add_pixel(const Ccoord_t<DimS> & /*pixel*/) {
     throw std::runtime_error
       ("this material needs pixels with a column vector of Euler angles");
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM, Int NbSlip>
   void MaterialCrystalPlasticityFinite<DimS, DimM, NbSlip>::
   add_pixel(const Ccoord_t<DimS> & pixel,
 	    const Eigen::Ref<Eigen::Array<Real, NbEuler, 1>> Euler){
     this->internal_fields.add_pixel(pixel);
     this->EulerField.push_back(Euler);
   }
 
 
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM, Int NbSlip>
+  void MaterialCrystalPlasticityFinite<DimS, DimM, NbSlip>::
+  initialise() {
+    if (not this->is_initialised) {
+      Parent::initialise();
+      using T2_t = Eigen::Matrix<Real, DimM, DimM>;
+      this->FpField.get_map().current() = T2_t::Identity();
+
+      using ColArray_t = Eigen::Matrix<Real, NbSlip, 1>;
+      this->GammaDotField.get_map().current() = ColArray_t::Zero();
+      this->TauYField.get_map().current() = ColArray_t::Constant(this->tau_y0);
+
+      this->GammaField.set_zero();
+
+      this->FpField.cycle();
+      this->GammaDotField.cycle();
+      this->TauYField.cycle();
+
+    }
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM, Int NbSlip>
+  void MaterialCrystalPlasticityFinite<DimS, DimM, NbSlip>::
+  save_history_variables() {
+    auto GammaMap = this->GammaField.get_map();
+    auto GammaDotMap = this->GammaDotField.get_map();
+
+    for (auto && tup: akantu::zip(GammaMap, GammaDotMap)) {
+      auto & gamma = std::get<0>(tup);
+      auto & gamma_dot = std::get<1>(tup);
+      gamma += .5 * this->delta_t * (gamma_dot.current() + gamma_dot.old());
+    }
+
+    this->FpField.cycle();
+    this->GammaDotField.cycle();
+    this->TauYField.cycle();
+
+  }
+
   /* ---------------------------------------------------------------------- */
   template class MaterialCrystalPlasticityFinite<  twoD,   twoD,  7>;
   template class MaterialCrystalPlasticityFinite<  twoD, threeD, 12>;
   template class MaterialCrystalPlasticityFinite<threeD, threeD, 12>;
 
 
 } // muSpectre
diff --git a/src/materials/material_crystal_plasticity_finite.hh b/src/materials/material_crystal_plasticity_finite.hh
index e9be273..c5ce051 100644
--- a/src/materials/material_crystal_plasticity_finite.hh
+++ b/src/materials/material_crystal_plasticity_finite.hh
@@ -1,454 +1,464 @@
 /**
  * @file   material_crystal_plasticity_finite.hh
  *
  * @author Till Junge <till.junge@altermail.ch>
  * @author Francesco Maresca <francesco.maresca@epfl.ch>
  *
  * @date   23 Feb 2018
  *
  * @brief finite strain crystal plasticity implementation
  *
  * Copyright © 2018 Till Junge, Francesco Maresca
  *
  * µ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 MATERIAL_CRYSTAL_PLASTICITY_FINITE_H
 #define MATERIAL_CRYSTAL_PLASTICITY_FINITE_H
 
 #include "materials/material_muSpectre_base.hh"
 #include "common/field.hh"
 #include "common/geometry.hh"
 #include "common/statefield.hh"
 #include "common/eigen_tools.hh"
 
 #include <Eigen/Dense>
 
 #include <cmath>
 
 namespace muSpectre {
 
   template <Dim_t DimS, Dim_t DimM, Int NbSlip>
   class MaterialCrystalPlasticityFinite;
 
   /**
    * traits for objective linear elasticity with eigenstrain
    */
   template <Dim_t DimS, Dim_t DimM, Int NbSlip>
   struct MaterialMuSpectre_traits<MaterialCrystalPlasticityFinite<DimS, DimM, NbSlip>>
   {
     //! global field collection
     using GFieldCollection_t = typename
       MaterialBase<DimS, DimM>::GFieldCollection_t;
 
     //! expected map type for strain fields
     using StrainMap_t = MatrixFieldMap<GFieldCollection_t, Real, DimM, DimM, true>;
     //! expected map type for stress fields
     using StressMap_t = MatrixFieldMap<GFieldCollection_t, Real, DimM, DimM>;
     //! expected map type for tangent stiffness fields
     using TangentMap_t = T4MatrixFieldMap<GFieldCollection_t, Real, DimM>;
 
     //! declare what type of strain measure your law takes as input
     constexpr static auto strain_measure{StrainMeasure::Gradient};
     //! declare what type of stress measure your law yields as output
     constexpr static auto stress_measure{StressMeasure::PK2};
 
     //! declare internal variables
     //! local field_collections used for internals
     using LFieldColl_t = LocalFieldCollection<DimS>;
     //! plastic deformation gradient Fₚ(t)
     using FpMap_t = StateFieldMap<MatrixFieldMap<LFieldColl_t, Real, DimM, DimM, false>>;
     //! plastic slip rates dγᵅ/dt
     using GammaDotMap_t = StateFieldMap<MatrixFieldMap<LFieldColl_t, Real, NbSlip, 1, false>>;
     //! critical resolved shear stresses (CRSS) τᵅy(t)
     using TauYMap_t = GammaDotMap_t;
-    //! plastic slips γᵅ(t)
-    using GammaMap_t = MatrixFieldMap<LFieldColl_t, Real, NbSlip, 1, false>;
     //! euler angles
     using EulerMap_t = ArrayFieldMap<LFieldColl_t, Real, (DimM==3) ? 3 : 1, 1, true>;
 
-    using InternalVariables = std::tuple<FpMap_t, GammaDotMap_t, TauYMap_t, GammaMap_t, EulerMap_t>;
+    using InternalVariables = std::tuple<FpMap_t, GammaDotMap_t, TauYMap_t, EulerMap_t>;
   };
 
   /**
    * implements finite strain crystal plasticity
    */
   template <Dim_t DimS, Dim_t DimM, Int NbSlip>
   class MaterialCrystalPlasticityFinite:
     public MaterialMuSpectre<MaterialCrystalPlasticityFinite<DimS, DimM, NbSlip>, DimS, DimM>
   {
   public:
     constexpr static Int NbEuler{(DimM==3) ? 3 : 1};
 
     //! base class
     using Parent = MaterialMuSpectre<MaterialCrystalPlasticityFinite, DimS, DimM>;
 
     //! type for stiffness tensor construction
     using Stiffness_t = Eigen::TensorFixedSize
       <Real, Eigen::Sizes<DimM, DimM, DimM, DimM>>;
 
     //! traits of this material
     using traits = MaterialMuSpectre_traits<MaterialCrystalPlasticityFinite>;
 
     //! Type of container used for storing eigenstrain
     using InternalVariables = typename traits::InternalVariables;
 
     //! Hooke's law implementation
     using Hooke = typename
       MatTB::Hooke<DimM,
                    typename traits::StrainMap_t::reference,
                    typename traits::TangentMap_t::reference>;
 
     //! Type in which plastic deformation gradient is referenced
     using Fp_ref = typename traits::FpMap_t::reference;
     //! Type in which slip rates are referenced
     using GammaDot_ref = typename traits::GammaDotMap_t::reference;
     //! Type in which CRSS are referenced
     using TauY_ref = typename traits::TauYMap_t::reference;
-    //! Type in which accumulated slips are referenced
-    using Gamma_ref = typename traits::GammaMap_t::reference;
     //! Type in which Euler angles are referenced
     using Euler_ref = typename traits::EulerMap_t::reference;
     //! Type in which slip directions and normals are given
     using SlipVecs = Eigen::Matrix<Real, NbSlip, DimM>;
     using SlipVecs_ref = Eigen::Ref<SlipVecs>;
     //! Type for second rank tensors
     using T2_t = Eigen::Matrix<Real, DimM, DimM>;
     //! Type for fourth rank tensors
     using T4_t = T4Mat<Real, DimM>;
 
     //! Default constructor
     MaterialCrystalPlasticityFinite() = delete;
 
     /**
      * Construct by name, Bulk modulus, Shear modulus, Reference slip
      * rate, Strain rate sensitivity, Initial CRSS, Initial hardening
      * modulus, CRSS saturation value, Hardening exponent,
      * Latent/Self-hardening ratio, Slip directions, Slip normals,
      * Plastic slip rate tolerance, Plastic slip rate loop maximum
      * iterations
      */
     MaterialCrystalPlasticityFinite(std::string name, Real bulk_m, Real shear_m,
                                     Real gamma_dot0, Real m_par, Real tau_y0,
                                     Real h0, Real s_infty, Real a_par,
                                     Real q_n, SlipVecs_ref Slip0,
                                     SlipVecs_ref Normal0, Real delta_t,
                                     Real tolerance=1.e-4, Int maxiter=20);
 
     //! Copy constructor
     MaterialCrystalPlasticityFinite(const MaterialCrystalPlasticityFinite &other) = delete;
 
     //! Move constructor
     MaterialCrystalPlasticityFinite(MaterialCrystalPlasticityFinite &&other) = delete;
 
     //! Destructor
     virtual ~MaterialCrystalPlasticityFinite() = default;
 
     //! Copy assignment operator
     MaterialCrystalPlasticityFinite& operator=(const MaterialCrystalPlasticityFinite &other) = delete;
 
     //! Move assignment operator
     MaterialCrystalPlasticityFinite& operator=(MaterialCrystalPlasticityFinite &&other) = delete;
 
     /**
      * evaluates second Piola-Kirchhoff stress given the Deformation Gradient
      */
     template <class s_t>
     inline decltype(auto) evaluate_stress(s_t && F, Fp_ref Fp,
                                           GammaDot_ref gamma_dot,
-                                          TauY_ref tau_y, Gamma_ref Gamma,
-                                          Euler_ref Euler);
+                                          TauY_ref tau_y, Euler_ref Euler);
 
     /**
      * evaluates both second Piola-Kirchhoff stress and stiffness given
      * the Deformation Gradient
      */
     template <class s_t>
     inline decltype(auto)
     evaluate_stress_tangent(s_t && F, Fp_ref Fp, GammaDot_ref gamma_dot,
-                            TauY_ref tau_y, Gamma_ref Gamma, Euler_ref Euler);
+                            TauY_ref tau_y, Euler_ref Euler);
 
     /**
      * return the internals tuple
      */
     InternalVariables & get_internals() {
       return this->internal_variables;};
 
     /**
      * overload add_pixel to write into internals
      */
     void add_pixel(const Ccoord_t<DimS> & pixel) override final;
 
     /**
      * overload add_pixel to write into eigenstrain
      */
     void add_pixel(const Ccoord_t<DimS> & pixel,
                    const Eigen::Ref<Eigen::Array<Real, NbEuler, 1>> Euler);
 
     /**
      * for introspection
      */
     constexpr static Dim_t get_NbSlip() {return NbSlip;}
 
+
+    /**
+     * set initial values for internal variables
+     */
+    void initialise() override;
+
+    /**
+     * material-specific update of internal (history) variables
+     */
+    void save_history_variables() override;
+
   protected:
     using LColl_t = typename traits::LFieldColl_t;
     //! Storage for F_p
     StateField<TensorField<LColl_t, Real, secondOrder, DimM>> FpField;
     //! Storage for dγ/dt
     StateField<MatrixField<LColl_t, Real, NbSlip, 1>> GammaDotField;
     //! Storage for τ_y
     StateField<MatrixField<LColl_t, Real, NbSlip, 1>> TauYField;
     //! Storage for γ
     MatrixField<LColl_t, Real, NbSlip, 1> & GammaField;
     //! Storage for Euler angles
     MatrixField<LColl_t, Real, NbEuler, 1> & EulerField;
     //! bulk modulus
     Real bulk_m;
     Real shear_m;
     Real gamma_dot0;
     Real m_par;
-    Real tau_y0; //// TODO: is not used? Francesco?
+    Real tau_y0;
     Real h0;
     Real s_infty;
     Real a_par;
     Real q_n;
     Real delta_t;
     Real tolerance;
     Int maxiter;
     //! Storage for slip directions
     alignas(16) Eigen::Matrix<Real, NbSlip, DimM> Slip0;
     //! Storage for slip plane normals
     alignas(16) Eigen::Matrix<Real, NbSlip, DimM> Normal0;
     T4_t C_el{};
 
     //! tuple for iterable internal variables
     InternalVariables internal_variables;
   private:
   };
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM, Int NbSlip>
   template <class s_t>
   decltype(auto)
   MaterialCrystalPlasticityFinite<DimS, DimM, NbSlip>::
   evaluate_stress(s_t && F, Fp_ref Fp, GammaDot_ref gamma_dot, TauY_ref tau_y,
-                  Gamma_ref Gamma, Euler_ref Euler) {
+                  Euler_ref Euler) {
     return std::get<0>(this->evaluate_stress_tangent(std::forward<s_t>(F),
-                                                     Fp, gamma_dot, tau_y, Gamma,
+                                                     Fp, gamma_dot, tau_y,
                                                      Euler));
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM, Int NbSlip>
   template <class s_t>
   decltype(auto)
   MaterialCrystalPlasticityFinite<DimS, DimM, NbSlip>::
   evaluate_stress_tangent(s_t && F, Fp_ref Fp, GammaDot_ref gamma_dot,
-                          TauY_ref tau_y, Gamma_ref /*Gamma*/,
-                          Euler_ref Euler) {
+                          TauY_ref tau_y, Euler_ref Euler) {
+    auto dot = [] (auto && a, auto && b) {return Matrices::dot<DimM>(a, b);};
+    auto ddot = [] (auto && a, auto && b) {return Matrices::ddot<DimM>(a, b);};
+
     Rotator<DimM> Rot(Euler);
     T2_t Floc{Rot.rotate(F)};
     std::array<T2_t, NbSlip> SchmidT;
     for (Int i{0}; i < NbSlip; ++i) {
       SchmidT[i] = this->Slip0.row(i).transpose() * this->Normal0.row(i);
     }
 
     // trial elastic deformation
     T2_t Fe_star{Floc*Fp.old().inverse()};
-    T2_t CGe_star{Fe_star.transpose()*Fe_star};
+    T2_t CGe_star{Fe_star.transpose()*Fe_star}; // elastic Cauchy-Green strain
     T2_t GLe_star{.5*(CGe_star - T2_t::Identity())};
     T2_t SPK_star{Matrices::tensmult(C_el,GLe_star)};
 
     using ColArray_t = Eigen::Array<Real, NbSlip, 1>;
     using ColMatrix_t = Eigen::Matrix<Real, NbSlip, 1>;
     ColArray_t tau_star;
-    // pl_corr is the plastic corrector
+    // pl_corr is the plastic corrector (Bᵅ)
     std::array<T2_t, NbSlip> pl_corr;
     using SlipMat_t = Eigen::Matrix<Real, NbSlip, NbSlip>;
     SlipMat_t pl_corr_proj;
 
     for (Int i{0}; i < NbSlip; ++i) {
       tau_star(i) = (CGe_star*SPK_star*SchmidT[i].transpose()).trace();
+
+      // eq (19)
       pl_corr[i] = Matrices::tensmult(C_el,.5*(CGe_star*SchmidT[i]+SchmidT[i].transpose()*CGe_star));
       for (Int j{0}; j < NbSlip; ++j) {
-        pl_corr_proj(i,j)=(Matrices::tensmult(C_el,pl_corr[i]*SchmidT[j].transpose())).trace();
+        pl_corr_proj(i,j)= ddot(CGe_star*pl_corr[i], SchmidT[j]);
       }
     }
 
     SlipMat_t I(SlipMat_t::Identity());
     auto && q_matrix{(I + this->q_n*(SlipMat_t::Ones() -I))};
 
     // residual on plastic slip rates
     gamma_dot.current() = gamma_dot.old();
     tau_y.current() = tau_y.old();
     ColArray_t tau_inc{tau_star};
 
     auto compute_gamma_dot = [this, &gamma_dot, &tau_inc, &tau_y] () {
       return gamma_dot.current().array()-this->gamma_dot0*(abs(tau_inc).array()/tau_y.current().array()).pow(this->m_par)*sign(tau_inc.array()); };
     ColArray_t res{compute_gamma_dot()};
 
     auto compute_h_matrix = [this, &q_matrix] (const ColMatrix_t & tau_y_temp) {
       auto && parens =
       (ColMatrix_t::Ones()-tau_y_temp/this->s_infty).array()
       .pow(this->a_par).matrix();
       return this->h0*parens.asDiagonal()*q_matrix; };
     ColArray_t s_dot_old{(compute_h_matrix(tau_y.current())*gamma_dot.old()).array()};
 
     SlipMat_t drdgammadot{SlipMat_t::Identity()};
     ColArray_t dr_stress{ColArray_t::Zero()};
 
     Int counter{};
 
     while (abs(res).maxCoeff()/this->gamma_dot0 > tolerance){
       if(counter ++ > this->maxiter){
         throw std::runtime_error("Max. number of iteration for plastic slip reached without convergence");
       }
 
       dr_stress = abs(tau_inc).pow((1-this->m_par)/this->m_par)*tau_y.current().array().pow(-1/this->m_par)*sign(tau_inc);
       ColArray_t dr_hard{abs(tau_inc).pow(1/this->m_par)*tau_y.current().array().pow(-1-1/this->m_par)*sign(tau_inc)};
       drdgammadot = I+0.5*this->delta_t*this->gamma_dot0/this->m_par*dr_stress.matrix().asDiagonal()*pl_corr_proj.transpose()
         +0.5*this->delta_t*this->gamma_dot0/this->m_par*dr_hard.matrix().asDiagonal()*compute_h_matrix(tau_y.current())*Eigen::sign(gamma_dot.current().array()).matrix().asDiagonal();
       gamma_dot.current() -= drdgammadot.inverse() * res.matrix();
-      // TODO: Check with Francesco whether the transposition of this guy is correct
-      tau_inc = tau_star - (0.5*this->delta_t*(gamma_dot.current()+gamma_dot.old()).transpose()*pl_corr_proj.transpose()).array().transpose();
+
+      tau_inc = tau_star -
+        (0.5*this->delta_t*(gamma_dot.current() + gamma_dot.old()).transpose() *
+         pl_corr_proj).array().transpose();
 
       Int counter_h{};
       ColMatrix_t tau_y_temp{};
       do {
         if(counter_h ++ > this->maxiter){
           throw std::runtime_error("Max. number of iteration for hardening reached without convergence");
         }
         tau_y_temp = tau_y.current();
         tau_y.current() = tau_y.old() + 0.5*this->delta_t*(s_dot_old.matrix() + compute_h_matrix(tau_y_temp)*gamma_dot.current());
       } while ((tau_y.current() - tau_y_temp).norm() > tolerance);
 
       res = compute_gamma_dot();
     }
 
     T2_t SPK{SPK_star};
     for (Int i{0}; i < NbSlip; ++i) {
       SPK -= .5*delta_t*(gamma_dot.current()(i)+gamma_dot.old()(i))*pl_corr[i];
     }
 
     T2_t Lp{T2_t::Zero()};
     for (Int i{0}; i < NbSlip; ++i) {
       Lp += .5*(gamma_dot.current()(i)+gamma_dot.old()(i))*SchmidT[i];
     }
 
     Fp.current() = (T2_t::Identity()+this->delta_t*Lp)*Fp.old();
 
     T2_t PK2 = Rot.rotate_back(Fp.current().inverse()*SPK*Fp.current().inverse().transpose());
 
       // Stiffness matrix calculation begins
 
       // A4: elastic trial consistent tangent
 
     auto IRT = Matrices::Itrns<DimM>();
     auto I4 = Matrices::Iiden<DimM>();
     auto odot = [] (auto && T4, auto && T2) {
       T4_t Return_value(T4_t::Zero());
       for (Int i = 0; i < DimM; ++i) {
         for (Int j = 0; j < DimM; ++j) {
           for (Int k = 0; k < DimM; ++k) {
             for (Int l = 0; l < DimM; ++l) {
               for (Int m = 0; m < DimM; ++m) {
                 get(Return_value,i,j,k,l) += get(T4,i,m,k,l)*T2(m,j);
               }
             }
           }
         }
       }
       return Return_value;
     };
 
-    auto dot = [] (auto && a, auto && b) {return Matrices::dot<DimM>(a, b);};
-    auto ddot = [] (auto && a, auto && b) {return Matrices::ddot<DimM>(a, b);};
     T4_t dAdF{odot(dot(Fp.old().inverse().transpose(),IRT),Fe_star)+odot(dot(Fe_star.transpose(),I4),Fp.old().inverse())};
     T4_t A4{.5*ddot(C_el,dAdF)};
 
     // E4: Tangent of the projector
 
     T4_t E4{T4_t::Zero()};
     for (Int i{0}; i < NbSlip; ++i) {
       T4_t dBprojdF{odot(dAdF,SchmidT[i]) + dot(SchmidT[i].transpose(),dAdF)};
       E4 -= .5*this->delta_t*(gamma_dot.current()(i)+gamma_dot.old()(i))*ddot(C_el,dBprojdF);
     }
 
     // G4: Tangent of slip rate
 
     // dgammadot/dtau
     SlipMat_t dgammadotdtau{-drdgammadot.inverse()*(-this->gamma_dot0/this->m_par*dr_stress.matrix().asDiagonal())};
 
     T4_t G4p1{T4_t::Zero()};
 
     for (Int k{0}; k < NbSlip; ++k) {
       for (Int mu{0}; mu < NbSlip; ++mu) {
 
         G4p1 += Matrices::outer(pl_corr[k],dgammadotdtau(k,mu)*SchmidT[mu]);
 
       }
     }
 
     T4_t G4p2{odot(dAdF,SPK)+dot(CGe_star,(A4+E4).eval())};
     T4_t G4_RHS{.5*delta_t*G4p1*G4p2};
 
     auto xdot = [] (auto && T4, auto && T2) {
       T4_t Return_value(T4_t::Zero());
       for (Int i = 0; i < DimM; ++i) {
         for (Int j = 0; j < DimM; ++j) {
           for (Int k = 0; k < DimM; ++k) {
             for (Int l = 0; l < DimM; ++l) {
               for (Int m = 0; m < DimM; ++m) {
                 get(Return_value,i,j,k,l) += get(T4,i,j,m,l)*T2(m,k);
               }
             }
           }
         }
       }
       return Return_value;
     };
 
     T4_t G4_LHS{I4-.5*delta_t*xdot(G4p1,CGe_star)};
 
     T4_t G4{G4_LHS.inverse()*G4_RHS};
 
     // Plastic contribution to geometric tangent
 
     T4_t F4p1{T4_t::Zero()};
     T4_t F4p2{T4_t::Zero()};
 
     for (Int k{0}; k < NbSlip; ++k) {
       for (Int mu{0}; mu < NbSlip; ++mu) {
 
         F4p1 += Matrices::outer(SchmidT[k],dgammadotdtau(k,mu)*SchmidT[mu]);
         F4p2 += Matrices::outer(SchmidT[k].transpose(),dgammadotdtau(k,mu)*SchmidT[mu]);
 
       }
     }
 
     T4_t F4L{-.5*delta_t*odot(dot(Fp.old(),(F4p1*(A4+E4+G4)).eval()),SPK*Fp.current().inverse().transpose())};
     T4_t F4R{-.5*delta_t*(dot((Fp.current().inverse()*SPK).eval(),odot((F4p2*(A4+E4+G4)).eval(),Fp.old().inverse().transpose())))};
 
     T4_t K4{dot(F,Rot.rotate_back(F4L + odot(dot(Fp.current().inverse(),(A4+E4+G4).eval()),Fp.current().inverse().transpose())  + F4R))};
 
     return std::make_tuple(std::move(PK2), std::move(K4));
 
   }
 
 
 }  // muSpectre
 
 #endif /* MATERIAL_CRYSTAL_PLASTICITY_FINITE_H */
diff --git a/src/materials/materials_toolbox.hh b/src/materials/materials_toolbox.hh
index f293d8c..e62a956 100644
--- a/src/materials/materials_toolbox.hh
+++ b/src/materials/materials_toolbox.hh
@@ -1,741 +1,740 @@
 /**
  * @file   materials_toolbox.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   02 Nov 2017
  *
  * @brief  collection of common continuum mechanics tools
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #ifndef MATERIALS_TOOLBOX_H
 #define MATERIALS_TOOLBOX_H
 
 #include "common/common.hh"
 #include "common/tensor_algebra.hh"
 #include "common/eigen_tools.hh"
 #include "common/T4_map_proxy.hh"
 
 #include <Eigen/Dense>
 #include <unsupported/Eigen/MatrixFunctions>
 
 #include <exception>
 #include <sstream>
 #include <iostream>
 #include <tuple>
 #include <type_traits>
 
 namespace muSpectre {
 
   namespace MatTB {
     /**
      * thrown when generic materials-related runtime errors occur
      * (mostly continuum mechanics problems)
      */
     class MaterialsToolboxError:public std::runtime_error{
     public:
       //! constructor
       explicit MaterialsToolboxError(const std::string& what)
         :std::runtime_error(what){}
       //! constructor
       explicit MaterialsToolboxError(const char * what)
         :std::runtime_error(what){}
     };
 
 
     /* ---------------------------------------------------------------------- */
     /**
      * Flag used to designate whether the material should compute both stress
      * and tangent moduli or only stress
      */
     enum class NeedTangent {
       yes, //!< compute both stress and tangent moduli
       no   //!< compute only stress
     };
 
     /**
      * struct used to determine the exact type of a tuple of references obtained
      * when a bunch of iterators over fiel_maps are dereferenced and their
      * results are concatenated into a tuple
      */
     template <class... T>
     struct ReferenceTuple {
       //! use this type
       using type = std::tuple<typename T::reference ...>;
     };
 
     /**
      * specialisation for tuples
      */
     //template <>
     template <class... T>
     struct ReferenceTuple<std::tuple<T...>> {
       //! use this type
       using type = typename ReferenceTuple<T...>::type;
     };
 
     /**
      * helper type for ReferenceTuple
      */
     template <class... T>
     using ReferenceTuple_t = typename ReferenceTuple<T...>::type;
 
     /* ---------------------------------------------------------------------- */
     namespace internal {
 
       /** Structure for functions returning one strain measure as a
        *  function of another
        **/
       template <StrainMeasure In, StrainMeasure Out = In>
       struct ConvertStrain {
         static_assert((In == StrainMeasure::Gradient) ||
                       (In == StrainMeasure::Infinitesimal),
                       "This situation makes me suspect that you are not using "
                       "MatTb as intended. Disable this assert only if you are "
                       "sure about what you are doing.");
 
         //! returns the converted strain
         template <class Strain_t>
         inline static decltype(auto)
         compute(Strain_t&& input) {
           // transparent case, in which no conversion is required:
           // just a perfect forwarding
           static_assert
             ((In == Out),
              "This particular strain conversion is not implemented");
           return std::forward<Strain_t>(input);
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /** Specialisation for getting Green-Lagrange strain from the
           transformation gradient
           E = ¹/₂ (C - I) = ¹/₂ (Fᵀ·F - I)
       **/
       template <>
       struct ConvertStrain<StrainMeasure::Gradient, StrainMeasure::GreenLagrange> {
 
         //! returns the converted strain
         template <class Strain_t>
         inline static decltype(auto)
         compute(Strain_t&& F) {
           return .5*(F.transpose()*F - Strain_t::PlainObject::Identity());
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /** Specialisation for getting Left Cauchy-Green strain from the
           transformation gradient
           B = F·Fᵀ = V²
       **/
       template <>
       struct ConvertStrain<StrainMeasure::Gradient, StrainMeasure::LCauchyGreen> {
 
         //! returns the converted strain
         template <class Strain_t>
         inline static decltype(auto)
         compute(Strain_t&& F) {
           return F*F.transpose();
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /** Specialisation for getting Right Cauchy-Green strain from the
           transformation gradient
           C = Fᵀ·F = U²
       **/
       template <>
       struct ConvertStrain<StrainMeasure::Gradient, StrainMeasure::RCauchyGreen> {
 
         //! returns the converted strain
         template <class Strain_t>
         inline static decltype(auto)
         compute(Strain_t&& F) {
           return F.transpose()*F;
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /** Specialisation for getting logarithmic (Hencky) strain from the
           transformation gradient
           E₀ = ¹/₂ ln C = ¹/₂ ln (Fᵀ·F)
       **/
       template <>
       struct ConvertStrain<StrainMeasure::Gradient, StrainMeasure::Log> {
 
         //! returns the converted strain
         template <class Strain_t>
         inline static decltype(auto)
         compute(Strain_t&& F) {
           constexpr Dim_t dim{EigenCheck::tensor_dim<Strain_t>::value};
           return (.5*logm(Eigen::Matrix<Real, dim, dim>{F.transpose()*F})).eval();
         }
       };
 
     }  // internal
 
     /* ---------------------------------------------------------------------- */
     //! set of functions returning one strain measure as a function of
     //! another
     template <StrainMeasure In, StrainMeasure Out,
               class Strain_t>
     decltype(auto) convert_strain(Strain_t && strain) {
       return internal::ConvertStrain<In, Out>::compute(std::move(strain));
     };
 
 
 
     /* ---------------------------------------------------------------------- */
     namespace internal {
 
       /** Structure for functions returning PK1 stress from other stress measures
        **/
       template <Dim_t Dim,
                 StressMeasure StressM,
                 StrainMeasure StrainM>
       struct PK1_stress {
 
         //! returns the converted stress
         template <class Strain_t, class Stress_t>
         inline static decltype(auto)
         compute(Strain_t && /*strain*/, Stress_t && /*stress*/) {
           // the following test always fails to generate a compile-time error
           static_assert((StressM == StressMeasure::Cauchy) &&
                         (StressM == StressMeasure::PK1),
                         "The requested Stress conversion is not implemented. "
                         "You either made a programming mistake or need to "
                         "implement it as a specialisation of this function. "
                         "See PK2stress<PK1,T1, T2> for an example.");
         }
 
         //! returns the converted stress and stiffness
         template <class Strain_t, class Stress_t, class Tangent_t>
         inline static decltype(auto)
         compute(Strain_t && /*strain*/, Stress_t && /*stress*/,
                 Tangent_t && /*stiffness*/) {
           // the following test always fails to generate a compile-time error
           static_assert((StressM == StressMeasure::Cauchy) &&
                         (StressM == StressMeasure::PK1),
                         "The requested Stress conversion is not implemented. "
                         "You either made a programming mistake or need to "
                         "implement it as a specialisation of this function. "
                         "See PK2stress<PK1,T1, T2> for an example.");
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /** Specialisation for the transparent case, where we already
           have PK1 stress
        **/
       template <Dim_t Dim, StrainMeasure StrainM>
       struct PK1_stress<Dim, StressMeasure::PK1, StrainM>:
         public PK1_stress<Dim, StressMeasure::no_stress_,
                           StrainMeasure::no_strain_> {
 
         //! returns the converted stress
         template <class Strain_t, class Stress_t>
         inline static decltype(auto)
         compute(Strain_t && /*dummy*/, Stress_t && P) {
           return std::forward<Stress_t>(P);
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /** Specialisation for the transparent case, where we already have PK1
           stress *and* stiffness is given with respect to the transformation
           gradient
        **/
       template <Dim_t Dim>
       struct PK1_stress<Dim, StressMeasure::PK1, StrainMeasure::Gradient>:
         public PK1_stress<Dim, StressMeasure::PK1,
                           StrainMeasure::no_strain_> {
 
         //! base class
         using Parent = PK1_stress<Dim, StressMeasure::PK1,
                                   StrainMeasure::no_strain_>;
         using Parent::compute;
 
         //! returns the converted stress and stiffness
         template <class Strain_t, class Stress_t, class Tangent_t>
         inline static decltype(auto)
         compute(Strain_t && /*dummy*/, Stress_t && P, Tangent_t && K) {
           return std::make_tuple(std::forward<Stress_t>(P),
                                  std::forward<Tangent_t>(K));
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /**
        * Specialisation for the case where we get material stress
        * (Piola-Kirchhoff-2, PK2)
        */
       template <Dim_t Dim, StrainMeasure StrainM>
       struct PK1_stress<Dim, StressMeasure::PK2, StrainM>:
         public PK1_stress<Dim, StressMeasure::no_stress_,
                           StrainMeasure::no_strain_> {
 
         //! returns the converted stress
         template <class Strain_t, class Stress_t>
         inline static decltype(auto)
         compute(Strain_t && F, Stress_t && S) {
           return F*S;
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /**
        * Specialisation for the case where we get material stress
        * (Piola-Kirchhoff-2, PK2) derived with respect to
        * Green-Lagrange strain
        */
       template <Dim_t Dim>
       struct PK1_stress<Dim, StressMeasure::PK2, StrainMeasure::GreenLagrange>:
         public PK1_stress<Dim, StressMeasure::PK2,
                           StrainMeasure::no_strain_> {
         //! base class
         using Parent = PK1_stress<Dim, StressMeasure::PK2,
                                   StrainMeasure::no_strain_>;
         using Parent::compute;
 
         //! returns the converted stress and stiffness
         template <class Strain_t, class Stress_t, class Tangent_t>
         inline static decltype(auto)
         compute(Strain_t && F, Stress_t && S, Tangent_t && C) {
           using T4 = typename std::remove_reference_t<Tangent_t>::PlainObject;
           using Tmap = T4MatMap<Real, Dim>;
           T4 K;
           Tmap Kmap{K.data()};
           K.setZero();
 
           for (int i = 0; i < Dim; ++i) {
             for (int m = 0; m < Dim; ++m) {
               for (int n = 0; n < Dim; ++n) {
                 get(Kmap,i,m,i,n) += S(m,n);
                 for (int j = 0; j < Dim; ++j) {
                   for (int r = 0; r < Dim; ++r) {
                     for (int s = 0; s < Dim; ++s) {
                       get(Kmap,i,m,j,n) += F(i,r)*get(C,r,m,n,s)*(F(j,s));
                     }
                   }
                 }
               }
             }
           }
           auto && P = compute(std::forward<Strain_t>(F),
                               std::forward<Stress_t>(S));
           return std::make_tuple(std::move(P), std::move(K));
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /**
        * Specialisation for the case where we get material stress
        * (Piola-Kirchhoff-2, PK2) derived with respect to
        * the placement Gradient (F)
        */
-      // TODO: check with Francesco
       template <Dim_t Dim>
       struct PK1_stress<Dim, StressMeasure::PK2, StrainMeasure::Gradient>:
         public PK1_stress<Dim, StressMeasure::PK2,
                           StrainMeasure::no_strain_> {
         //! base class
         using Parent = PK1_stress<Dim, StressMeasure::PK2,
                                   StrainMeasure::no_strain_>;
         using Parent::compute;
 
         //! returns the converted stress and stiffness
         template <class Strain_t, class Stress_t, class Tangent_t>
         inline static decltype(auto)
         compute(Strain_t && F, Stress_t && S, Tangent_t && C) {
           using T4 = typename std::remove_reference_t<Tangent_t>::PlainObject;
           using Tmap = T4MatMap<Real, Dim>;
           T4 K;
           Tmap Kmap{K.data()};
           K.setZero();
 
           for (int i = 0; i < Dim; ++i) {
             for (int m = 0; m < Dim; ++m) {
               for (int n = 0; n < Dim; ++n) {
                 get(Kmap,i,m,i,n) += S(m,n);
                 for (int j = 0; j < Dim; ++j) {
                   for (int r = 0; r < Dim; ++r) {
                     for (int s = 0; s < Dim; ++s) {
                       get(Kmap,i,m,j,n) += F(i,r)*get(C,r,m,j,n);
                     }
                   }
                 }
               }
             }
           }
           auto && P = compute(std::forward<Strain_t>(F),
                               std::forward<Stress_t>(S));
           return std::make_tuple(std::move(P), std::move(K));
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /**
        * Specialisation for the case where we get Kirchhoff stress (τ)
        */
       template <Dim_t Dim, StrainMeasure StrainM>
       struct PK1_stress<Dim, StressMeasure::Kirchhoff, StrainM>:
         public PK1_stress<Dim, StressMeasure::no_stress_,
                           StrainMeasure::no_strain_> {
 
         //! returns the converted stress
         template <class Strain_t, class Stress_t>
         inline static decltype(auto)
         compute(Strain_t && F, Stress_t && tau) {
           return tau*F.inverse().transpose();
         }
       };
 
       /* ---------------------------------------------------------------------- */
       /**
        * Specialisation for the case where we get Kirchhoff stress (τ) derived
        * with respect to Gradient
        */
       template <Dim_t Dim>
       struct PK1_stress<Dim, StressMeasure::Kirchhoff, StrainMeasure::Gradient>:
         public PK1_stress<Dim, StressMeasure::Kirchhoff,
                           StrainMeasure::no_strain_> {
 
         //! short-hand
         using Parent = PK1_stress<Dim, StressMeasure::Kirchhoff,
                                   StrainMeasure::no_strain_>;
         using Parent::compute;
 
         //! returns the converted stress and stiffness
         template <class Strain_t, class Stress_t, class Tangent_t>
         inline static decltype(auto)
         compute(Strain_t && F, Stress_t && tau, Tangent_t && C) {
           using T4 = typename std::remove_reference_t<Tangent_t>::PlainObject;
           using Tmap = T4MatMap<Real, Dim>;
           T4 K;
           Tmap Kmap{K.data()};
           K.setZero();
           auto && F_inv{F.inverse()};
           for (int i = 0; i < Dim; ++i) {
             for (int m = 0; m < Dim; ++m) {
               for (int n = 0; n < Dim; ++n) {
                 for (int j = 0; j < Dim; ++j) {
                   for (int r = 0; r < Dim; ++r) {
                     for (int s = 0; s < Dim; ++s) {
                       get(Kmap,i,m,j,n) += F_inv(i,r)*get(C,r,m,n,s);
                     }
                   }
                 }
               }
             }
           }
           auto && P = tau * F_inv.transpose();
           return std::make_tuple(std::move(P), std::move(K));
         }
       };
 
     }  // internal
 
     /* ---------------------------------------------------------------------- */
     //! set of functions returning an expression for PK2 stress based on
     template <StressMeasure StressM, StrainMeasure StrainM,
               class Stress_t, class Strain_t>
     decltype(auto) PK1_stress(Strain_t && strain, Stress_t && stress) {
       constexpr Dim_t dim{EigenCheck::tensor_dim<Strain_t>::value};
       static_assert((dim == EigenCheck::tensor_dim<Stress_t>::value),
                     "Stress and strain tensors have differing dimensions");
       return internal::PK1_stress<dim, StressM, StrainM>::compute
         (std::forward<Strain_t>(strain),
          std::forward<Stress_t>(stress));
     };
 
     /* ---------------------------------------------------------------------- */
     //! set of functions returning an expression for PK2 stress based on
     template <StressMeasure StressM, StrainMeasure StrainM,
               class Stress_t, class Strain_t, class Tangent_t>
     decltype(auto) PK1_stress(Strain_t  && strain,
                               Stress_t  && stress,
                               Tangent_t && tangent) {
       constexpr Dim_t dim{EigenCheck::tensor_dim<Strain_t>::value};
       static_assert((dim == EigenCheck::tensor_dim<Stress_t>::value),
                     "Stress and strain tensors have differing dimensions");
       static_assert((dim == EigenCheck::tensor_4_dim<Tangent_t>::value),
                     "Stress and tangent tensors have differing dimensions");
 
       return internal::PK1_stress<dim, StressM, StrainM>::compute
         (std::forward<Strain_t>(strain),
          std::forward<Stress_t>(stress),
          std::forward<Tangent_t>(tangent));
     };
 
     namespace internal {
 
       //! Base template for elastic modulus conversion
       template <ElasticModulus Out, ElasticModulus In1, ElasticModulus In2>
       struct Converter {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real& /*in1*/,
                                       const Real& /*in2*/) {
           // if no return has happened until now, the conversion is not
           // implemented yet
           static_assert((In1 == In2),
                         "This conversion has not been implemented yet, please add "
                         "it here below as a specialisation of this function "
                         "template. Check "
                         "https://en.wikipedia.org/wiki/Lam%C3%A9_parameters for "
                         "the formula.");
           return 0;
         }
       };
 
       /**
        * Spectialisation for when the output is the first input
        */
       template <ElasticModulus Out, ElasticModulus In>
       struct Converter<Out, Out, In> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & A, const Real & /*B*/) {
           return A;
         }
       };
 
       /**
        * Spectialisation for when the output is the second input
        */
       template <ElasticModulus Out, ElasticModulus In>
       struct Converter<Out, In, Out> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & /*A*/, const Real & B) {
           return B;
         }
       };
 
       /**
        * Specialisation μ(E, ν)
        */
       template <>
       struct Converter<ElasticModulus::Shear,
                        ElasticModulus::Young,
                        ElasticModulus::Poisson> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & E, const Real & nu) {
           return E/(2*(1+nu));
         }
       };
 
       /**
        * Specialisation λ(E, ν)
        */
       template <>
       struct Converter<ElasticModulus::lambda,
                        ElasticModulus::Young,
                        ElasticModulus::Poisson> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & E, const Real & nu) {
           return E*nu/((1+nu)*(1-2*nu));
         }
       };
 
       /**
        * Specialisation K(E, ν)
        */
       template <>
       struct Converter<ElasticModulus::Bulk,
                        ElasticModulus::Young,
                        ElasticModulus::Poisson> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & E, const Real & nu) {
           return E/(3*(1-2*nu));
         }
       };
 
       /**
        * Specialisation E(K, µ)
        */
       template <>
       struct Converter<ElasticModulus::Young,
                        ElasticModulus::Bulk,
                        ElasticModulus::Shear> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & K, const Real & G) {
         return 9*K*G/(3*K+G);
         }
       };
 
       /**
        * Specialisation ν(K, µ)
        */
       template <>
       struct Converter<ElasticModulus::Poisson,
                        ElasticModulus::Bulk,
                        ElasticModulus::Shear> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & K, const Real & G) {
           return (3*K - 2*G) /(2*(3*K + G));
         }
       };
 
       /**
        * Specialisation E(λ, µ)
        */
       template <>
       struct Converter<ElasticModulus::Young,
                        ElasticModulus::lambda,
                        ElasticModulus::Shear> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & lambda, const Real & G) {
           return G * (3*lambda + 2*G)/(lambda + G);
         }
       };
 
       /**
        * Specialisation λ(K, µ)
        */
       template <>
       struct Converter<ElasticModulus::lambda,
                        ElasticModulus::Bulk,
                        ElasticModulus::Shear> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & K, const Real & mu) {
           return K - 2. * mu / 3.;
         }
       };
 
     }  // internal
 
     /**
      * allows the conversion from any two distinct input moduli to a
      * chosen output modulus
      */
     template <ElasticModulus Out, ElasticModulus In1, ElasticModulus In2>
     inline constexpr Real convert_elastic_modulus(const Real& in1,
                                                   const Real& in2) {
       // enforcing sanity
       static_assert((In1 != In2),
                     "The input modulus types cannot be identical");
 
       // enforcing independence from order in which moduli are supplied
       constexpr bool inverted{In1 > In2};
       using Converter = std::conditional_t<inverted,
                                            internal::Converter<Out, In2, In1>,
                                            internal::Converter<Out, In1, In2>>;
       if (inverted) {
         return Converter::compute(std::move(in2),
                                   std::move(in1));
       } else {
         return Converter::compute(std::move(in1),
                                   std::move(in2));
       }
     }
 
 
     //! static inline implementation of Hooke's law
     template <Dim_t Dim, class Strain_t, class Tangent_t>
     struct Hooke {
       /**
        * compute Lamé's first constant
        * @param young: Young's modulus
        * @param poisson: Poisson's ratio
        */
       inline static constexpr Real
       compute_lambda(const Real & young, const Real & poisson) {
         return convert_elastic_modulus<ElasticModulus::lambda,
                                        ElasticModulus::Young,
                                        ElasticModulus::Poisson>(young, poisson);
       }
 
       /**
        * compute Lamé's second constant (i.e., shear modulus)
        * @param young: Young's modulus
        * @param poisson: Poisson's ratio
        */
       inline static constexpr Real
       compute_mu(const Real & young, const Real & poisson) {
         return convert_elastic_modulus<ElasticModulus::Shear,
                                        ElasticModulus::Young,
                                        ElasticModulus::Poisson>(young, poisson);
       }
 
       /**
        * compute the bulk modulus
        * @param young: Young's modulus
        * @param poisson: Poisson's ratio
        */
       inline static constexpr Real
       compute_K(const Real & young, const Real & poisson) {
         return convert_elastic_modulus<ElasticModulus::Bulk,
                                        ElasticModulus::Young,
                                        ElasticModulus::Poisson>(young, poisson);
       }
 
       /**
        * compute the stiffness tensor
        * @param lambda: Lamé's first constant
        * @param mu: Lamé's second constant (i.e., shear modulus)
        */
       inline static Eigen::TensorFixedSize<Real, Eigen::Sizes<Dim, Dim, Dim, Dim>>
       compute_C(const Real & lambda, const Real & mu) {
         return lambda*Tensors::outer<Dim>(Tensors::I2<Dim>(),Tensors::I2<Dim>()) +
           2*mu*Tensors::I4S<Dim>();
       }
 
       /**
        * compute the stiffness tensor
        * @param lambda: Lamé's first constant
        * @param mu: Lamé's second constant (i.e., shear modulus)
        */
       inline static T4Mat<Real, Dim>
       compute_C_T4(const Real & lambda, const Real & mu) {
         return lambda*Matrices::Itrac<Dim>() + 2*mu*Matrices::Isymm<Dim>();
       }
 
       /**
        * return stress
        * @param lambda: First Lamé's constant
        * @param mu: Second Lamé's constant (i.e. shear modulus)
        * @param E: Green-Lagrange or small strain tensor
        */
       template <class s_t>
       inline static decltype(auto)
       evaluate_stress(const Real & lambda, const Real & mu, s_t && E) {
         return E.trace()*lambda * Strain_t::Identity() + 2*mu*E;
       }
 
       /**
        * return stress and tangent stiffness
        * @param lambda: First Lamé's constant
        * @param mu: Second Lamé's constant (i.e. shear modulus)
        * @param E: Green-Lagrange or small strain tensor
        * @param C: stiffness tensor (Piola-Kirchhoff 2 (or σ) w.r.t to `E`)
        */
       template <class s_t>
       inline static decltype(auto)
       evaluate_stress(const Real & lambda, const Real & mu,
                       Tangent_t && C, s_t && E) {
         return std::make_tuple
           (std::move(evaluate_stress(lambda, mu, std::move(E))),
            std::move(C));
       }
     };
 
   }  // MatTB
 
 }  // muSpectre
 
 #endif /* MATERIALS_TOOLBOX_H */