diff --git a/doc/dev-docs/source/Doxyfile b/doc/dev-docs/source/Doxyfile
index 441cc25..aa66c89 100644
--- a/doc/dev-docs/source/Doxyfile
+++ b/doc/dev-docs/source/Doxyfile
@@ -1,24 +1,26 @@
 PROJECT_NAME           = muSpectre
 ## this is for read the docs, which builds in-place
 INPUT                  = ../../../src
 ## this overwrites the INPUT path for local builds
 @INCLUDE               = input_def
 RECURSIVE              = YES
 
+EXTRACT_ALL            = YES
+
 SOURCE_BROWSER         = YES
 
 GENERATE_HTML          = YES
 GENERATE_LATEX         = NO
 GENERATE_XML           = YES
 ## this is for read the docs, which builds in-place
 XML_OUTPUT             = doxygenxml
 ## this overwrites the XML_OUTPUT path for local builds
 @INCLUDE               =  xml_output_def
 XML_PROGRAMLISTING     = YES
 
 ALIASES                = "rst=\verbatim embed:rst"
 ALIASES               += "endrst=\endverbatim"
 
 QUIET                  = YES
 WARNINGS               = YES
 WARN_IF_UNDOCUMENTED   = YES
diff --git a/doc/dev-docs/source/Tutorials.rst b/doc/dev-docs/source/GettingStarted.rst
similarity index 100%
copy from doc/dev-docs/source/Tutorials.rst
copy to doc/dev-docs/source/GettingStarted.rst
diff --git a/doc/dev-docs/source/NewMaterial.rst b/doc/dev-docs/source/NewMaterial.rst
new file mode 100644
index 0000000..de42d64
--- /dev/null
+++ b/doc/dev-docs/source/NewMaterial.rst
@@ -0,0 +1,191 @@
+Writing a New Constitutive Law
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The abstraction for a constitutive law in **µSpectre** is the ``Material``, and new such materials must inherit from the class :cpp:class:`muSpectre::MaterialBase`. Most often, however, it will be most convenient to inherit from the derived class :cpp:class:`muSpectre::MaterialMuSpectre`, as it implements a lot of the machinery that is commonly used by constitutive laws. This section describes how to implement a new constitutive law with internal variables (sometimes also called state variables). The example material implemented here is :cpp:class:`MaterialTutorial`, an objective linear elastic law with a distribution of eigenstrains as internal variables. The constitutive law is defined by the relationship between material (or second Piola-Kirchhoff) stress :math:`\mathbf{S}` and Green-Lagrange strain :math:`\mathbf{E}`
+
+.. math::
+   :nowrap:
+
+   \begin{align}
+      \mathbf{S} &= \mathbb C:\left(\mathbf{E}-\mathbf{e}\right), \\
+      S_{ij} &= C_{ijkl}\left(E_{kl}-e_{kl}\right), \\
+   \end{align}
+
+where :math:`\mathbb C` is the elastic stiffness tensor and :math:`\mathbf e` is the local eigenstrain. Note that the implementation sketched out here is the most convenient to quickly get started with using **µSpectre**, but not the most efficient one. For a most efficient implementation, refer to the implementation of :cpp:class:`muSpectre::MaterialLinearElastic2`.
+
+The :cpp:class:`muSpectre::MaterialMuSpectre` class
+***************************************************
+
+The class :cpp:class:`muSpectre::MaterialMuSpectre` is defined in ``material_muSpectre_base.hh`` and takes three template parameters;
+
+#. ``class Material`` is a `CRTP <https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern>`_ and names the material inheriting from it. The reason for this construction is that we want to avoid virtual method calls from :cpp:class:`muSpectre::MaterialMuSpectre` to its derived classes. Rather, we want :cpp:class:`muSpectre::MaterialMuSpectre` to be able to call methods of the inheriting class directly without runtime overhead.
+#. ``Dim_t DimS`` defines the number of spatial dimensions of the problem, i.e., whether we are dealing with a two- or three-dimensional grid of pixels/voxels.
+#. ``Dim_t DimM`` defines the number of dimensions of our material description. This value will typically be the same as ``DimS``, but in cases like generalised plane strain, we can for instance have a three three-dimensional material response in a two-dimensional pixel grid.
+
+The main job of :cpp:class:`muSpectre::MaterialMuSpectre` is to
+
+#. loop over all pixels to which this material has been assigned, transform the global gradient :math:`\mathbf{F}` (or small strain tensor :math:`\boldsymbol\varepsilon`) into the new material's required strain measure (e.g., the Green-Lagrange strain tensor :math:`\mathbf{E}`),
+#. for each pixel evaluate the constitutive law by calling its ``evaluate_stress`` (computes the stress response) or ``evaluate_stress_tangent`` (computes both stress and consistent tangent) method with the local strain and internal variables, and finally
+#. transform the stress (and possibly tangent) response from the material's stress measure into first Piola-Kirchhoff stress :math:`\mathbf{P}` (or Cauchy stress :math:`\boldsymbol\sigma` in small strain).
+
+In order to use these facilities, the new material needs to inherit from :cpp:class:`muSpectre::MaterialMuSpectre` (where we calculation of the response) and specialise the type :cpp:class:`muSpectre::MaterialMuSpectre_traits` (where we tell :cpp:class:`muSpectre::MaterialMuSpectre` how to use the new material). These two steps are described here for our example material.
+
+Specialising the :cpp:class:`muSpectre::MaterialMuSpectre_traits` structure
+***************************************************************************
+This structure is templated by the new material (in this case :cpp:class:`MaterialTutorial`) and needs to specify
+
+#. the types used to communicate per-pixel strains, stresses and stiffness tensors to the material (i.e., whether you want to get maps to `Eigen` matrices or raw pointers, or ...). Here we will use the convenient :cpp:type:`muSpectre::MatrixFieldMap` for strains and stresses, and :cpp:type:`muSpectre::T4MatrixFieldMap` for the stiffness. Look through the classes deriving from :cpp:type:`muSpectre::FieldMap` for all available options.
+#. the strain measure that is expected (e.g., gradient, Green-Lagrange strain, left Cauchy-Green strain, etc.). Here we will use Green-Lagrange strain. The options are defined by the enum :cpp:enum:`muSpectre::StrainMeasure`.
+#. the stress measure that is computed by the law (e.g., Cauchy, first Piola-Kirchhoff, etc,). Here, it will be first Piola-Kirchhoff stress. The available options are defined by the enum :cpp:enum:`muSpectre::StressMeasure`.
+
+Our traits look like this (assuming we are in the namespace ``muSpectre``::
+
+  template <Dim_t DimS, Dim_t DimM>
+  struct MaterialMuSpectre_traits<MaterialTutorial<DimS, DimM>>
+  {
+    //! global field collection
+    using GFieldCollection_t = typename
+      GlobalFieldCollection<DimS, DimM>;
+
+    //! 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::GreenLagrange};
+    //! declare what type of stress measure your law yields as output
+    constexpr static auto stress_measure{StressMeasure::PK2};
+
+    //! local field_collections used for internals
+    using LFieldColl_t = LocalFieldCollection<DimS, DimM>;
+    //! local strain type
+    using LStrainMap_t = MatrixFieldMap<LFieldColl_t, Real, DimM, DimM, true>;
+    //! elasticity with eigenstrain
+    using InternalVariables = std::tuple<LStrainMap_t>;
+
+  };
+
+Implementing the new material
+*****************************
+
+The new law needs to implement the methods ``add_pixel``, ``get_internals``, ``evaluate_stress``, and ``evaluate_stress_tangent``. Below is a commented example header::
+
+  template <Dim_t DimS, Dim_t DimM>
+  class MaterialTutorial:
+    public MaterialMuSpectre<MaterialTutorial<DimS, DimM>, DimS, DimM>
+  {
+  public:
+    //! traits of this material
+    using traits = MaterialMuSpectre_traits<MaterialTutorial>;
+
+    //! Type of container used for storing eigenstrain
+    using InternalVariables = typename traits::InternalVariables;
+
+    //! Construct by name, Young's modulus and Poisson's ratio
+    MaterialTutorial(std::string name, Real young, Real poisson);
+
+    /**
+     * evaluates second Piola-Kirchhoff stress given the Green-Lagrange
+     * strain (or Cauchy stress if called with a small strain tensor)
+     */
+    template <class s_t, class eigen_s_t>
+    inline decltype(auto) evaluate_stress(s_t && E, eigen_s_t && E_eig);
+
+    /**
+     * evaluates both second Piola-Kirchhoff stress and stiffness given
+     * the Green-Lagrange strain (or Cauchy stress and stiffness if
+     * called with a small strain tensor)
+     */
+    template <class s_t, class eigen_s_t>
+    inline decltype(auto)
+    evaluate_stress_tangent(s_t &&  E, eigen_s_t && E_eig);
+
+    /**
+     * return the internals tuple (needed by `muSpectre::MaterialMuSpectre`)
+     */
+    InternalVariables & get_internals() {
+      return this->internal_variables;};
+
+    /**
+     * overload add_pixel to write into eigenstrain
+     */
+    void add_pixel(const Ccoord_t<DimS> & pixel,
+                   const Eigen::Matrix<Real, DimM, DimM> & E_eig);
+
+  protected:
+    //! stiffness tensor
+    T4Mat<Real, DimM> C;
+    //! storage for eigenstrain
+    using Field_t =
+      TensorField<LocalFieldCollection<DimS,DimM>, Real, secondOrder, DimM>;
+    Field_t & eigen_field; //!< field of eigenstrains
+    //! tuple for iterable eigen_field
+    InternalVariables internal_variables;
+  private:
+  };
+
+A possible implementation for the constructor would be::
+
+  template <Dim_t DimS, Dim_t DimM>
+  MaterialTutorial<DimS, DimM>::MaterialTutorial(std::string name,
+                                                 Real young,
+                                                 Real poisson)
+    :MaterialMuSpectre<MaterialTutorial, DimS, DimM>(name) {
+
+    // Lamé parameters
+    Real lambda{young*poisson/((1+poisson)*(1-2*poisson))};
+    Real mu{young/(2*(1+poisson))};
+
+    // Kronecker delta
+    Eigen::Matrix<Real, DimM, DimM> del{Eigen::Matrix<Real, DimM, DimM>::Identity()};
+
+
+    // fill the stiffness tensor
+    this->C.setZero();
+    for (Dim_t i = 0; i < DimM; ++i) {
+      for (Dim_t j = 0; j < DimM; ++j) {
+        for (Dim_t k = 0; k < DimM; ++k) {
+          for (Dim_t l = 0; l < DimM; ++l) {
+            get(this->C, i, j, k, l) += (lambda * del(i,j)*del(k,l) +
+                                         mu * (del(i,k)*del(j,l) + del(i,l)*del(j,k)));
+          }
+        }
+      }
+    }
+  }
+
+as an exercise, you could check how :cpp:class:`muSpectre::MaterialLinearElastic1` uses **µSpectre**'s materials toolbox (in namespace ``MatTB``) to compute :math:`\mathbb C` in a much more convenient fashion. The evaluation of the stress could be (here, we make use of the ``Matrices`` namespace that defines common tensor algebra operations)::
+
+  template <Dim_t DimS, Dim_t DimM>
+  template <class s_t, class eigen_s_t>
+  decltype(auto)
+  MaterialTutorial<DimS, DimM>::
+  evaluate_stress(s_t && E, eigen_s_t && E_eig) {
+    return Matrices::tens_mult(this->C, E-E_eig);
+  }
+
+
+
+The remaining two methods are straight-forward::
+
+  template <Dim_t DimS, Dim_t DimM>
+  template <class s_t, class eigen_s_t>
+  decltype(auto)
+  MaterialTutorial<DimS, DimM>::
+  evaluate_stress_tangent(s_t && E, eigen_s_t && E_eig) {
+    return return std::make_tuple
+          (evaluate_stress(E, E_eig),
+           this->C);
+  }
+
+  template <Dim_t DimS, Dim_t DimM>
+  InternalVariables &
+  MaterialTutorial<DimS, DimM>::get_internals() {
+    return this->internal_variables;
+  }
+
+
+Note that the methods ``evaluate_stress`` and ``evaluate_stress_tangent`` need to be in the header, as both their input parameter types and output type depend on the compile-time context.
diff --git a/doc/dev-docs/source/Reference.rst b/doc/dev-docs/source/Reference.rst
index 1c092a6..2352b07 100644
--- a/doc/dev-docs/source/Reference.rst
+++ b/doc/dev-docs/source/Reference.rst
@@ -1,13 +1,7 @@
 .. _reference:
 
 Reference
 ---------
 
 .. doxygenindex::
    :project: µSpectre
-
-.. Doxygennamespace
-   ----------------
-.. .. doxygennamespace:: muSpectre
-      :project: µSpectre
-      :outline:
diff --git a/doc/dev-docs/source/Tutorials.rst b/doc/dev-docs/source/Tutorials.rst
index 9606557..a59140e 100644
--- a/doc/dev-docs/source/Tutorials.rst
+++ b/doc/dev-docs/source/Tutorials.rst
@@ -1,92 +1,6 @@
-Getting Started
-~~~~~~~~~~~~~~~
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents:
 
-Obtaining **µSpectre**
-**********************
-
-**µSpectre** is hosted on a git repository on `c4science`_. To clone it, run
-
-.. code-block:: sh
-
-   $ git clone https://c4science.ch/source/muSpectre.git
-
-or if you prefer identifying yourself using a public ssh-key, run
-
-.. code-block:: bash
-
-   $ git clone ssh://git@c4science.ch/source/muSpectre.git
-
-The latter option requires you to have a user account on c4science (`create <https://c4science.ch/auth/start/>`_).
-
-
-.. _c4science: https://c4science.ch
-
-
-Building **µSpectre**
-*********************
-
-You can compile **µSpectre** using  `CMake <https://cmake.org/>`_. The current (and possibly incomplete list of) dependencies are
-
-- `Boost unit test framework <http://www.boost.org/doc/libs/1_66_0/libs/test/doc/html/index.html>`_
-- `FFTW <http://www.fftw.org>`_
-- `Sphinx <http://www.sphinx-doc.org>`_ and `Breathe <https://breathe.readthedocs.io>`_
-
-**µSpectre** requires a relatively modern compiler as it makes heavy use of C++14 features. It has successfully been compiled and tested using the following compilers under Linux
-
-- gcc-7.2
-- gcc-6.4
-- gcc-5.4
-- clang-5.0
-- clang-4.0
-
-and using clang-4.0 under MacOS.
-
-It does *not* compile on Intel's most recent compiler, as it is still lacking some C++14 support. Work-arounds are planned, but will have to wait for someone to pick up the `task <https://c4science.ch/T1852>`_.
-
-To compile, create a build folder and configure the CMake project. If you do this in the folder you cloned in the previous step, it can look for instance like this:
-
-.. code-block:: sh
-
-   $ mkdir build-release
-   $ cd build-release
-   $ ccmake ..
-
-Then, set the build type to ``Release`` to produce optimised code. µSpectre makes heavy use of expression templates, so optimisation is paramount. (As an example, the performance difference between code compiled in ``Debug`` and ``Release`` is about a factor 40 in simple linear elasticity.)
-
-Finally, compile the library and the tests by running
-
-.. code-block:: sh
-
-   $ make -j <NB-OF-PROCESSES>
-
-.. warning::
-
-   When using the ``-j`` option to compile, be aware that compiling µSpectre uses quite a bit of RAM. If your machine start swapping at compile time, reduce the number of parallel compilations
-
-
-Running **µSpectre**
-********************
-
-The easiest and intended way of using **µSpectre** is through its Python bindings. The following simple example computes the response of a two-dimensional stretched periodic RVE cell. The cell consist of a soft matrix with a circular hard inclusion.
-
-.. literalinclude:: ../../../bin/tutorial_example.py
-   :language: python
-
-More examples both both python and c++ executables can be found in the ``/bin`` folder.
-
-Getting help
-************
-
-µSpectre is in a very early stage of development and the documentation is currently spotty. Also, there is no FAQ page yet. If you run into trouble, please contact us on the `µSpectre chat room <https://c4science.ch/Z69>`_ and someone will answer as soon as possible. You can also check the API :ref:`reference`.
-
-
-Reporting Bugs
-**************
-
-If you think you found a bug, you are probably right. Please report it! The preferred way is for you to create a task on `µSpectre's workboard <https://c4science.ch/project/board/1447/>`_ and assign it to user ``junge``. Include steps to reproduce the bug if possible. Someone will answer as soon as possible.
-
-
-Contribute
-**********
-
-We welcome contributions both for new features and bug fixes. New features must be documented and have unit tests. Please submit contributions for review as `Arcanist revisions <https://secure.phabricator.com/book/phabricator/article/arcanist/>`_. More detailed guidelines for submissions will follow soonᵀᴹ.
+   GettingStarted
+   NewMaterial
diff --git a/doc/dev-docs/source/index.rst b/doc/dev-docs/source/index.rst
index 204bf2f..164c633 100644
--- a/doc/dev-docs/source/index.rst
+++ b/doc/dev-docs/source/index.rst
@@ -1,39 +1,40 @@
 
 
 .. image:: ../logo_flat.png
 
 µSpectre, FFT-based Homogenisation without Linear Reference Medium
 ==================================================================
 
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents:
+
+   Summary
+   Tutorials
+   Reference
+   License
+
 µSpectre is free software; you can redistribute it and/or
 modify it under the :term:``s 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.
 
-.. toctree::
-   :maxdepth: 2
-   :caption: Contents:
-
-   Summary
-   Tutorials
-   Reference
-   License
 
 
 
 Indices and tables
 ==================
 
 * :ref:`genindex`
 * :ref:`modindex`
 * :ref:`search`
diff --git a/src/common/field_map_matrixlike.hh b/src/common/field_map_matrixlike.hh
index 332bc3f..75d7a41 100644
--- a/src/common/field_map_matrixlike.hh
+++ b/src/common/field_map_matrixlike.hh
@@ -1,402 +1,402 @@
 /**
  * @file   field_map_matrixlike.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   26 Sep 2017
  *
  * @brief  Eigen-Matrix and -Array maps over strongly typed fields
  *
  * 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 FIELD_MAP_MATRIXLIKE_H
 #define FIELD_MAP_MATRIXLIKE_H
 
 #include "common/field_map_base.hh"
 #include "common/T4_map_proxy.hh"
 
 #include <Eigen/Dense>
 
 #include <memory>
 
 namespace muSpectre {
 
   namespace internal {
 
     /* ---------------------------------------------------------------------- */
     /**
      * lists all matrix-like types consideres by
      * `muSpectre::internal::MatrixLikeFieldMap`
      */
     enum class Map_t{
       Matrix,  //!< for wrapping `Eigen::Matrix`
       Array,   //!< for wrapping `Eigen::Array`
       T4Matrix //!< for wrapping `Eigen::T4Matrix`
     };
 
     /**
      * traits structure to define the name shown when a
      * `muSpectre::MatrixLikeFieldMap` output into an ostream
      */
     template<Map_t map_type>
     struct NameWrapper {
     };
 
     /// specialisation for `muSpectre::ArrayFieldMap`
     template<>
     struct NameWrapper<Map_t::Array> {
       //! string to use for printing
       static std::string field_info_root() {return "Array";}
     };
 
     /// specialisation for `muSpectre::MatrixFieldMap`
     template<>
     struct NameWrapper<Map_t::Matrix> {
       //! string to use for printing
       static std::string field_info_root() {return "Matrix";}
     };
 
     /// specialisation for `muSpectre::T4MatrixFieldMap`
     template<>
     struct NameWrapper<Map_t::T4Matrix> {
       //! string to use for printing
       static std::string field_info_root() {return "T4Matrix";}
     };
 
     /* ---------------------------------------------------------------------- */
     /*!
      * base class for maps of matrices, arrays and fourth-order
      * tensors mapped onty matrices
      *
      * It should never be necessary to call directly any of the
      * constructors if this class, but rather use the template aliases
      * `muSpectre::ArrayFieldMap`, `muSpectre::MatrixFieldMap`, and
      * `muSpectre::T4MatrixFieldMap`
      */
     template <class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type, bool ConstField>
     class MatrixLikeFieldMap: public FieldMap<FieldCollection,
                                               typename EigenArray::Scalar,
                                               EigenArray::SizeAtCompileTime,
                                               ConstField>
     {
     public:
       //! base class
       using Parent = FieldMap<FieldCollection,
                               typename EigenArray::Scalar,
                               EigenArray::SizeAtCompileTime, ConstField>;
       using T_t = EigenPlain; //!< plain Eigen type to map
       //! cell coordinates type
       using Ccoord = Ccoord_t<FieldCollection::spatial_dim()>;
       using value_type = EigenArray; //!< stl conformance
       using const_reference = EigenConstArray; //!< stl conformance
       //! stl conformance
       using reference = std::conditional_t<ConstField,
                                            const_reference,
                                            value_type>; // since it's a resource handle
       using size_type = typename Parent::size_type; //!< stl conformance
       using pointer = std::unique_ptr<EigenArray>; //!< stl conformance
       using TypedField = typename Parent::TypedField; //!< stl conformance
       using Field = typename TypedField::Base; //!< stl conformance
       //! stl conformance
       using const_iterator= typename Parent::template iterator<MatrixLikeFieldMap, true>;
       //! stl conformance
       using iterator = std::conditional_t<
         ConstField,
         const_iterator,
         typename Parent::template iterator<MatrixLikeFieldMap>>;
       using reverse_iterator = std::reverse_iterator<iterator>; //!< stl conformance
        //! stl conformance
       using const_reverse_iterator = std::reverse_iterator<const_iterator>;
       //! stl conformance
       friend iterator;
 
       //! Default constructor
       MatrixLikeFieldMap() = delete;
 
       /**
        * Constructor using a (non-typed) field. Compatibility is enforced at
        * runtime.  This should not be a performance concern, as this constructor
        * will not be called in anny inner loops (if used as intended).
        */
       template <bool isntConst=!ConstField>
       MatrixLikeFieldMap(std::enable_if_t<isntConst, Field &> field);
       /**
        * Constructor using a (non-typed) field. Compatibility is enforced at
        * runtime.  This should not be a performance concern, as this constructor
        * will not be called in anny inner loops (if used as intended).
        */
       template <bool isConst=ConstField>
       MatrixLikeFieldMap(std::enable_if_t<isConst, const Field &> field);
 
       /**
        * Constructor using a typed field. Compatibility is enforced
        * statically. It is not always possible to call this constructor, as the
        * type of the field might not be known at compile time.
        */
       template<class FC, typename T2, Dim_t NbC>
       MatrixLikeFieldMap(TypedSizedFieldBase<FC, T2, NbC> & field);
 
       //! Copy constructor
       MatrixLikeFieldMap(const MatrixLikeFieldMap &other) = default;
 
       //! Move constructor
       MatrixLikeFieldMap(MatrixLikeFieldMap &&other) = default;
 
       //! Destructor
       virtual ~MatrixLikeFieldMap()  = default;
 
       //! Copy assignment operator
       MatrixLikeFieldMap& operator=(const MatrixLikeFieldMap &other) = delete;
 
       //! Move assignment operator
       MatrixLikeFieldMap& operator=(MatrixLikeFieldMap &&other) = delete;
 
       //! Assign a matrixlike value to every entry
       template <class Derived>
       inline MatrixLikeFieldMap & operator=(const Eigen::EigenBase<Derived> & val);
 
       //! give human-readable field map type
       inline std::string info_string() const override final;
 
       //! member access
       inline reference operator[](size_type index);
       //! member access
       inline reference operator[](const Ccoord& ccoord);
 
       //! member access
       inline const_reference operator[](size_type index) const;
       //! member access
       inline const_reference operator[](const Ccoord& ccoord) const;
 
       //! return an iterator to head of field for ranges
       inline iterator begin(){return iterator(*this);}
       //! return an iterator to head of field for ranges
       inline const_iterator cbegin() const {return const_iterator(*this);}
       //! return an iterator to head of field for ranges
       inline const_iterator begin() const {return this->cbegin();}
       //! return an iterator to tail of field for ranges
       inline iterator end(){return iterator(*this, false);};
       //! return an iterator to tail of field for ranges
       inline const_iterator cend() const {return const_iterator(*this, false);}
       //! return an iterator to tail of field for ranges
       inline const_iterator end() const {return this->cend();}
 
       //! evaluate the average of the field
       inline T_t mean() const;
 
     protected:
       //! for sad, legacy iterator use
       inline pointer ptr_to_val_t(size_type index);
       const static std::string field_info_root; //!< for printing and debug
     private:
     };
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type,  bool ConstField>
     const std::string MatrixLikeFieldMap<FieldCollection, EigenArray,
                                          EigenConstArray, EigenPlain, map_type,
                                          ConstField>::
     field_info_root{NameWrapper<map_type>::field_info_root()};
 
     /* ---------------------------------------------------------------------- */
     template<class FieldCollection, class EigenArray, class EigenConstArray,
              class EigenPlain, Map_t map_type,  bool ConstField>
     template <bool isntConst>
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
                        map_type, ConstField>::
     MatrixLikeFieldMap(std::enable_if_t<isntConst, Field &> field)
       :Parent(field) {
       static_assert((isntConst != ConstField),
                     "let the compiler deduce isntConst, this is a SFINAE "
                     "parameter");
       this->check_compatibility();
     }
 
     /* ---------------------------------------------------------------------- */
     template<class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type,  bool ConstField>
     template <bool isConst>
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
                        map_type, ConstField>::
     MatrixLikeFieldMap(std::enable_if_t<isConst, const Field &> field)
       :Parent(field) {
       static_assert((isConst == ConstField),
                     "let the compiler deduce isntConst, this is a SFINAE "
                     "parameter");
       this->check_compatibility();
     }
 
     /* ---------------------------------------------------------------------- */
     template<class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type,  bool ConstField>
     template<class FC, typename T2, Dim_t NbC>
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                        EigenPlain, map_type, ConstField>::
     MatrixLikeFieldMap(TypedSizedFieldBase<FC, T2, NbC> & field)
       :Parent(field) {
     }
 
     /* ---------------------------------------------------------------------- */
     //! human-readable field map type
     template<class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type,  bool ConstField>
     std::string
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                        EigenPlain, map_type, ConstField>::
     info_string() const {
       std::stringstream info;
       info << field_info_root << "("
            << typeid(typename EigenArray::value_type).name() << ", "
            << EigenArray::RowsAtCompileTime << "x"
            << EigenArray::ColsAtCompileTime << ")";
       return info.str();
     }
 
     /* ---------------------------------------------------------------------- */
     //! member access
     template<class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type,  bool ConstField>
     typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                                 EigenPlain, map_type, ConstField>::reference
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                        EigenPlain, map_type, ConstField>::
     operator[](size_type index) {
       return reference(this->get_ptr_to_entry(index));
     }
 
     template<class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type,  bool ConstField>
     typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                                 EigenPlain, map_type, ConstField>::reference
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
                        map_type, ConstField>::
     operator[](const Ccoord & ccoord) {
       size_t index{};
       index = this->collection.get_index(ccoord);
       return reference(this->get_ptr_to_entry(std::move(index)));
     }
 
     /* ---------------------------------------------------------------------- */
     //! member access
     template<class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type,  bool ConstField>
     typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                                 EigenPlain, map_type, ConstField>::const_reference
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
                        map_type, ConstField>::
     operator[](size_type index) const {
       return const_reference(this->get_ptr_to_entry(index));
     }
 
     template<class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type,  bool ConstField>
     typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                                 EigenPlain, map_type, ConstField>::const_reference
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
                        map_type, ConstField>::
     operator[](const Ccoord & ccoord) const{
       size_t index{};
       index = this->collection.get_index(ccoord);
       return const_reference(this->get_ptr_to_entry(std::move(index)));
     }
 
     //----------------------------------------------------------------------------//
     template <class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type,  bool ConstField>
     template <class Derived>
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                        EigenPlain, map_type, ConstField> &
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                        EigenPlain, map_type, ConstField>::
     operator=(const Eigen::EigenBase<Derived> & val) {
       for (auto && entry: *this) {
         entry = val;
       }
       return *this;
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type, bool ConstField>
     typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                                 EigenPlain, map_type, ConstField>::T_t
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                        EigenPlain, map_type, ConstField>::
     mean() const {
       T_t mean{T_t::Zero()};
       for (auto && val: *this) {
         mean += val;
       }
       mean *= 1./Real(this->size());
       return mean;
     }
 
     /* ---------------------------------------------------------------------- */
     template<class FieldCollection, class EigenArray, class EigenConstArray,
               class EigenPlain, Map_t map_type,  bool ConstField>
     typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
                                 EigenPlain, map_type, ConstField>::pointer
     MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
                        map_type, ConstField>::
     ptr_to_val_t(size_type index) {
       return std::make_unique<value_type>
         (this->get_ptr_to_entry(std::move(index)));
     }
 
   }  // internal
 
-  /* ---------------------------------------------------------------------- */
-  //! short-hand for an Eigen matrix map as iterate
-  template <class FieldCollection, typename T, Dim_t NbRows, Dim_t NbCols,
-            bool ConstField=false>
-  using MatrixFieldMap = internal::MatrixLikeFieldMap
-    <FieldCollection,
-     Eigen::Map<Eigen::Matrix<T, NbRows, NbCols>>,
-     Eigen::Map<const Eigen::Matrix<T, NbRows, NbCols>>,
-     Eigen::Matrix<T, NbRows, NbCols>,
-     internal::Map_t::Matrix, ConstField>;
-
   /* ---------------------------------------------------------------------- */
   //! short-hand for an Eigen matrix map as iterate
   template <class FieldCollection, typename T, Dim_t Dim,
             bool MapConst=false>
   using T4MatrixFieldMap = internal::MatrixLikeFieldMap
     <FieldCollection,
      T4MatMap<T, Dim, false>,
      T4MatMap<T, Dim, true>,
      T4Mat<T, Dim>,
      internal::Map_t::T4Matrix, MapConst>;
 
+  /* ---------------------------------------------------------------------- */
+  //! short-hand for an Eigen matrix map as iterate
+  template <class FieldCollection, typename T, Dim_t NbRows, Dim_t NbCols,
+            bool ConstField=false>
+  using MatrixFieldMap = internal::MatrixLikeFieldMap
+    <FieldCollection,
+     Eigen::Map<Eigen::Matrix<T, NbRows, NbCols>>,
+     Eigen::Map<const Eigen::Matrix<T, NbRows, NbCols>>,
+     Eigen::Matrix<T, NbRows, NbCols>,
+     internal::Map_t::Matrix, ConstField>;
+
   /* ---------------------------------------------------------------------- */
   //! short-hand for an Eigen array map as iterate
   template <class FieldCollection, typename T, Dim_t NbRows, Dim_t NbCols,
             bool ConstField=false>
   using ArrayFieldMap = internal::MatrixLikeFieldMap
     <FieldCollection,
      Eigen::Map<Eigen::Array<T, NbRows, NbCols>>,
      Eigen::Map<const Eigen::Array<T, NbRows, NbCols>>,
      Eigen::Array<T, NbRows, NbCols>,
      internal::Map_t::Array, ConstField>;
 
 
 }  // muSpectre
 
 #endif /* FIELD_MAP_MATRIXLIKE_H */
diff --git a/src/common/field_map_tensor.hh b/src/common/field_map_tensor.hh
index cab4d16..0112757 100644
--- a/src/common/field_map_tensor.hh
+++ b/src/common/field_map_tensor.hh
@@ -1,287 +1,287 @@
 /**
  * @file   field_map_tensor.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   26 Sep 2017
  *
  * @brief  Defines an Eigen-Tensor map over strongly typed fields
  *
  * 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 FIELD_MAP_TENSOR_H
 #define FIELD_MAP_TENSOR_H
 
 #include "common/eigen_tools.hh"
 #include "common/field_map_base.hh"
 
 #include <sstream>
 #include <memory>
 
 namespace muSpectre {
 
 
   /* ---------------------------------------------------------------------- */
   /**
    * maps onto a `muSpectre::internal::TypedSizedFieldBase` and lets
    * you iterate over it in the form ot
    * `Eigen::TensorMap<TensorFixedSize<...>>`
    */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
             bool ConstField=false>
   class TensorFieldMap: public
     internal::FieldMap<FieldCollection,
                        T,
                        SizesByOrder<order, dim>::Sizes::total_size,
                        ConstField>
   {
   public:
     //! base class
     using Parent = internal::FieldMap<FieldCollection, T,
                                       TensorFieldMap::nb_components, ConstField>;
     //! cell coordinates type
     using Ccoord = Ccoord_t<FieldCollection::spatial_dim()>;
     //! static tensor size
     using Sizes = typename SizesByOrder<order, dim>::Sizes;
     //! plain Eigen type
     using T_t = Eigen::TensorFixedSize<T, Sizes>;
     using value_type = Eigen::TensorMap<T_t>; //!< stl conformance
     using const_reference = Eigen::TensorMap<const T_t>; //!< stl conformance
     //! stl conformance
     using reference = std::conditional_t<ConstField,
                                          const_reference,
                                          value_type>; // since it's a resource handle
     using size_type = typename Parent::size_type; //!< stl conformance
     using pointer = std::unique_ptr<value_type>; //!< stl conformance
     using TypedField = typename Parent::TypedField; //!< field to map
     //! polymorphic base field type (for debug and python)
     using Field = typename TypedField::Base;
     //! stl conformance
     using const_iterator = typename Parent::template iterator<TensorFieldMap, true>;
     //! stl conformance
     using iterator = std::conditional_t<
       ConstField,
       const_iterator,
       typename Parent::template iterator<TensorFieldMap>>;
     //! stl conformance
     using reverse_iterator = std::reverse_iterator<iterator>;
     //! stl conformance
     using const_reverse_iterator = std::reverse_iterator<const_iterator>;
     //! stl conformance
     friend iterator;
 
     //! Default constructor
     TensorFieldMap() = delete;
 
     //! constructor
     template <bool isntConst=!ConstField>
     TensorFieldMap(std::enable_if_t<isntConst, Field &> field);
     //! constructor
     template <bool isConst=ConstField>
     TensorFieldMap(std::enable_if_t<isConst, const Field &> field);
 
     //! Copy constructor
     TensorFieldMap(const TensorFieldMap &other) = default;
 
     //! Move constructor
     TensorFieldMap(TensorFieldMap &&other) = default;
 
     //! Destructor
     virtual ~TensorFieldMap() = default;
 
     //! Copy assignment operator
     TensorFieldMap& operator=(const TensorFieldMap &other) = delete;
 
     //! Assign a matrixlike value to every entry
     inline TensorFieldMap & operator=(const T_t & val);
 
     //! Move assignment operator
     TensorFieldMap& operator=(TensorFieldMap &&other) = delete;
 
     //! give human-readable field map type
     inline std::string info_string() const override final;
 
     //! member access
     inline reference operator[](size_type index);
     //! member access
     inline reference operator[](const Ccoord & ccoord);
 
     //! member access
     inline const_reference operator[](size_type index) const;
     //! member access
     inline const_reference operator[](const Ccoord& ccoord) const;
 
     //! return an iterator to first entry of field
     inline iterator begin(){return iterator(*this);}
     //! return an iterator to first entry of field
     inline const_iterator cbegin() const {return const_iterator(*this);}
     //! return an iterator to first entry of field
     inline const_iterator begin() const {return this->cbegin();}
     //! return an iterator past the last entry of field
     inline iterator end(){return iterator(*this, false);}
     //! return an iterator past the last entry of field
     inline const_iterator cend() const {return const_iterator(*this, false);}
     //! return an iterator past the last entry of field
     inline const_iterator end() const {return this->cend();}
 
     //! evaluate the average of the field
     inline T_t mean() const;
   protected:
     //! for sad, legacy iterator use
     inline pointer ptr_to_val_t(size_type index);
   private:
   };
 
   /* ---------------------------------------------------------------------- */
   template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
            bool ConstField>
   template <bool isntConst>
   TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   TensorFieldMap(std::enable_if_t<isntConst, Field &> field)
     :Parent(field) {
     static_assert((isntConst != ConstField),
                   "let the compiler deduce isntConst, this is a SFINAE "
                   "parameter");
     this->check_compatibility();
   }
 
   /* ---------------------------------------------------------------------- */
   template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
            bool ConstField>
   template <bool isConst>
   TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   TensorFieldMap(std::enable_if_t<isConst, const Field &> field)
     :Parent(field) {
     static_assert((isConst == ConstField),
-                  "let the compiler deduce isntConst, this is a SFINAE "
+                  "let the compiler deduce isConst, this is a SFINAE "
                   "parameter");
     this->check_compatibility();
   }
 
   /* ---------------------------------------------------------------------- */
   //! human-readable field map type
   template<class FieldCollection, typename T, Dim_t order, Dim_t dim, bool ConstField>
   std::string
   TensorFieldMap<FieldCollection, T, order, dim, ConstField>::info_string() const {
     std::stringstream info;
     info << "Tensor(" << typeid(T).name() << ", " << order
          << "_o, " << dim << "_d)";
     return info.str();
   }
 
   /* ---------------------------------------------------------------------- */
   //! member access
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
             bool ConstField>
   typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::reference
   TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   operator[](size_type index) {
     auto && lambda = [this, &index](auto&&...tens_sizes) {
       return reference(this->get_ptr_to_entry(index), tens_sizes...);
     };
     return call_sizes<order, dim>(lambda);
   }
 
   template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
            bool ConstField>
   typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::reference
   TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   operator[](const Ccoord & ccoord) {
     auto && index = this->collection.get_index(ccoord);
     auto && lambda = [this, &index](auto&&...sizes) {
       return reference(this->get_ptr_to_entry(index), sizes...);
     };
     return call_sizes<order, dim>(lambda);
   }
 
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
             bool ConstField>
   typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   const_reference
   TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   operator[](size_type index) const {
     // Warning: due to a inconsistency in Eigen's API, tensor maps
     // cannot be constructed from a const ptr, hence this nasty const
     // cast :(
     auto && lambda = [this, &index](auto&&...tens_sizes) {
       return const_reference(const_cast<T*>(this->get_ptr_to_entry(index)),
                    tens_sizes...);
     };
     return call_sizes<order, dim>(lambda);
   }
 
   template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
            bool ConstField>
   typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   const_reference
   TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   operator[](const Ccoord & ccoord) const {
     auto && index = this->collection.get_index(ccoord);
     auto && lambda = [this, &index](auto&&...sizes) {
       return const_reference(const_cast<T*>(this->get_ptr_to_entry(index)),
                              sizes...);
     };
     return call_sizes<order, dim>(lambda);
   }
 
   /* ---------------------------------------------------------------------- */
   template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
            bool ConstField>
   TensorFieldMap<FieldCollection, T, order, dim, ConstField> &
   TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   operator=(const T_t & val) {
     for (auto && tens: *this) {
       tens = val;
     }
     return *this;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
             bool ConstField>
   typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::T_t
   TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   mean() const {
     T_t mean{T_t::Zero()};
     for (auto && val: *this) {
       mean += val;
     }
     mean *= 1./Real(this->size());
     return mean;
   }
 
   /* ---------------------------------------------------------------------- */
   //! for sad, legacy iterator use. Don't use unless you have to.
   template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
            bool ConstField>
   typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::pointer
   TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
   ptr_to_val_t(size_type index) {
     auto && lambda = [this, &index](auto&&... tens_sizes) {
       return std::make_unique<value_type>(this->get_ptr_to_entry(index),
                                        tens_sizes...);
     };
     return call_sizes<order, dim>(lambda);
   }
 
 
 }  // muSpectre
 
 
 #endif /* FIELD_MAP_TENSOR_H */
diff --git a/src/materials/material_linear_elastic2.hh b/src/materials/material_linear_elastic2.hh
index 5575437..74c7ee3 100644
--- a/src/materials/material_linear_elastic2.hh
+++ b/src/materials/material_linear_elastic2.hh
@@ -1,205 +1,199 @@
 /**
  * @file   material_linear_elastic2.hh
  *
  * @author Till Junge <till.junge@altermail.ch>
  *
  * @date   03 Feb 2018
  *
  * @brief linear elastic material with imposed eigenstrain and its
  *        type traits. Uses the MaterialMuSpectre facilities to keep it
  *        simple
  *
  * 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 MATERIAL_LINEAR_ELASTIC_EIGENSTRAIN_H
 #define MATERIAL_LINEAR_ELASTIC_EIGENSTRAIN_H
 
 #include "materials/material_linear_elastic1.hh"
 #include "common/field.hh"
 
 #include <Eigen/Dense>
 
 namespace muSpectre {
 
   template <Dim_t DimS, Dim_t DimM>
   class MaterialLinearElastic2;
 
   /**
    * traits for objective linear elasticity with eigenstrain
    */
   template <Dim_t DimS, Dim_t DimM>
   struct MaterialMuSpectre_traits<MaterialLinearElastic2<DimS, DimM>>
   {
     //! 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::GreenLagrange};
     //! declare what type of stress measure your law yields as output
     constexpr static auto stress_measure{StressMeasure::PK2};
 
     //! local field_collections used for internals
     using LFieldColl_t = LocalFieldCollection<DimS, DimM>;
     //! local strain type
     using LStrainMap_t = MatrixFieldMap<LFieldColl_t, Real, DimM, DimM, true>;
     //! elasticity with eigenstrain
     using InternalVariables = std::tuple<LStrainMap_t>;
 
   };
 
   /**
    * implements objective linear elasticity with an eigenstrain per pixel
    */
   template <Dim_t DimS, Dim_t DimM>
   class MaterialLinearElastic2:
     public MaterialMuSpectre<MaterialLinearElastic2<DimS, DimM>, DimS, DimM>
   {
   public:
 
     //! base class
     using Parent = MaterialMuSpectre<MaterialLinearElastic2, DimS, DimM>;
-    /**
-     * type used to determine whether the
-     * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only
-     * stresses or also tangent stiffnesses
-     */
-    using NeedTangent = typename Parent::NeedTangent;
 
     //! 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<MaterialLinearElastic2>;
 
     //! 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>;
 
     //! reference to any type that casts to a matrix
     using StrainTensor = Eigen::Ref<Eigen::Matrix<Real, DimM, DimM>>;
     //! Default constructor
     MaterialLinearElastic2() = delete;
 
     //! Construct by name, Young's modulus and Poisson's ratio
     MaterialLinearElastic2(std::string name, Real young, Real poisson);
 
 
     //! Copy constructor
     MaterialLinearElastic2(const MaterialLinearElastic2 &other) = delete;
 
     //! Move constructor
     MaterialLinearElastic2(MaterialLinearElastic2 &&other) = delete;
 
     //! Destructor
     virtual ~MaterialLinearElastic2() = default;
 
     //! Copy assignment operator
     MaterialLinearElastic2& operator=(const MaterialLinearElastic2 &other) = delete;
 
     //! Move assignment operator
     MaterialLinearElastic2& operator=(MaterialLinearElastic2 &&other) = delete;
 
     /**
      * evaluates second Piola-Kirchhoff stress given the Green-Lagrange
      * strain (or Cauchy stress if called with a small strain tensor)
      */
     template <class s_t, class eigen_s_t>
     inline decltype(auto) evaluate_stress(s_t && E, eigen_s_t && E_eig);
 
     /**
      * evaluates both second Piola-Kirchhoff stress and stiffness given
      * the Green-Lagrange strain (or Cauchy stress and stiffness if
      * called with a small strain tensor)
      */
     template <class s_t, class eigen_s_t>
     inline decltype(auto)
     evaluate_stress_tangent(s_t &&  E, eigen_s_t && E_eig);
 
     /**
      * return the internals tuple
      */
     InternalVariables & get_internals() {
       return this->internal_variables;};
 
     /**
      * overload add_pixel to write into eigenstrain
      */
     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 StrainTensor & E_eig);
 
   protected:
     //! linear material without eigenstrain used to compute response
     MaterialLinearElastic1<DimS, DimM> material;
     //! storage for eigenstrain
     using Field_t =
       TensorField<LocalFieldCollection<DimS,DimM>, Real, secondOrder, DimM>;
     Field_t & eigen_field; //!< field of eigenstrains
     //! tuple for iterable eigen_field
     InternalVariables internal_variables;
   private:
   };
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   template <class s_t, class eigen_s_t>
   decltype(auto)
   MaterialLinearElastic2<DimS, DimM>::
   evaluate_stress(s_t && E, eigen_s_t && E_eig) {
     return this->material.evaluate_stress(E-E_eig);
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   template <class s_t, class eigen_s_t>
   decltype(auto)
   MaterialLinearElastic2<DimS, DimM>::
   evaluate_stress_tangent(s_t && E, eigen_s_t && E_eig) {
     // using mat = Eigen::Matrix<Real, DimM, DimM>;
     // mat ecopy{E};
     // mat eig_copy{E_eig};
     // mat ediff{ecopy-eig_copy};
     // std::cout << "eidff - (E-E_eig)" << std::endl << ediff-(E-E_eig) << std::endl;
     // std::cout << "P1 <internal>" << std::endl << mat{std::get<0>(this->material.evaluate_stress_tangent(E-E_eig))} << "</internal>" << std::endl;
     // std::cout << "P2" << std::endl << mat{std::get<0>(this->material.evaluate_stress_tangent(std::move(ediff)))} << std::endl;
     return this->material.evaluate_stress_tangent(E-E_eig);
   }
 
 
 }  // muSpectre
 
 #endif /* MATERIAL_LINEAR_ELASTIC_EIGENSTRAIN_H */
diff --git a/src/solver/solver_cg.cc b/src/solver/solver_cg.cc
index 4692b78..4107b5d 100644
--- a/src/solver/solver_cg.cc
+++ b/src/solver/solver_cg.cc
@@ -1,128 +1,128 @@
 /**
  * @file   solver_cg.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   20 Dec 2017
  *
  * @brief  Implementation of cg solver
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "solver/solver_cg.hh"
 #include "solver/solver_error.hh"
 
 #include <iomanip>
 #include <cmath>
 #include <sstream>
 
 namespace muSpectre {
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   SolverCG<DimS, DimM>::SolverCG(Sys_t& sys, Real tol, Uint maxiter,
                                  bool verbose)
     :Parent(sys, tol, maxiter, verbose),
      r_k{make_field<Field_t>("residual r_k", this->collection)},
      p_k{make_field<Field_t>("search direction r_k", this->collection)},
      Ap_k{make_field<Field_t>("Effect of tangent A*p_k", this->collection)}
   {}
 
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void SolverCG<DimS, DimM>::solve(const Field_t & rhs,
                                    Field_t & x_f) {
     x_f.eigenvec() = this-> solve(rhs.eigenvec(), x_f.eigenvec());
   };
 
   //----------------------------------------------------------------------------//
   template <Dim_t DimS, Dim_t DimM>
   typename SolverCG<DimS, DimM>::SolvVectorOut
   SolverCG<DimS, DimM>::solve(const SolvVectorInC rhs, SolvVectorIn x_0) {
     // Following implementation of algorithm 5.2 in Nocedal's Numerical Optimization (p. 112)
 
     auto r = this->r_k.eigen();
     auto p = this->p_k.eigen();
     auto Ap = this->Ap_k.eigen();
-    auto x = typename Field_t::EigenMap(x_0.data(), r.rows(), r.cols());
+    typename Field_t::EigenMap x(x_0.data(), r.rows(), r.cols());
 
     // initialisation of algo
     r = this->sys.directional_stiffness_with_copy(x);
 
     r -= typename Field_t::ConstEigenMap(rhs.data(), r.rows(), r.cols());
     p = -r;
 
     this->converged = false;
     Real rdr = (r*r).sum();
     Real rhs_norm2 = rhs.squaredNorm();
     Real tol2 = ipow(this->tol,2)*rhs_norm2;
 
     size_t count_width{}; // for output formatting in verbose case
     if (this->verbose) {
       count_width = size_t(std::log10(this->maxiter))+1;
     }
 
     for (Uint i = 0;
          i < this->maxiter && (rdr > tol2 || i == 0);
          ++i, ++this->counter) {
       Ap = this->sys.directional_stiffness_with_copy(p);
 
       Real alpha = rdr/(p*Ap).sum();
 
       x += alpha * p;
       r += alpha * Ap;
 
       Real new_rdr = (r*r).sum();
       Real beta = new_rdr/rdr;
       rdr = new_rdr;
 
       if (this->verbose) {
         std::cout << "  at CG step " << std::setw(count_width) << i
                   << ": |r|/|b| = " << std::setw(15) << std::sqrt(rdr/rhs_norm2)
                   << ", cg_tol = " << this->tol << std::endl;
       }
 
       p = -r+beta*p;
     }
     if (rdr < tol2) {
       this->converged=true;
     } else {
       std::stringstream err {};
       err << " After " << this->counter << " steps, the solver "
           << " FAILED with  |r|/|b| = "
           << std::setw(15) << std::sqrt(rdr/rhs_norm2)
           << ", cg_tol = " << this->tol << std::endl;
       throw ConvergenceError("Conjugate gradient has not converged." + err.str());
     }
     return x_0;
   }
 
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   typename SolverCG<DimS, DimM>::Tg_req_t
   SolverCG<DimS, DimM>::get_tangent_req() const {
     return tangent_requirement;
   }
 
   template class SolverCG<twoD, twoD>;
   //template class SolverCG<twoD, threeD>;
   template class SolverCG<threeD, threeD>;
 }  // muSpectre