diff --git a/src/common/T4_map_proxy.hh b/src/common/T4_map_proxy.hh
index d151f66..3ea93c8 100644
--- a/src/common/T4_map_proxy.hh
+++ b/src/common/T4_map_proxy.hh
@@ -1,192 +1,193 @@
 /**
  * file   T4_map_proxy.hh
  *
  * @author Till Junge <till.junge@altermail.ch>
  *
  * @date   19 Nov 2017
  *
  * @brief  Map type to allow fourth-order tensor-like maps on 2D matrices
  *
  * @section LICENCE
  *
  * Copyright (C) 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include <type_traits>
 
 #include <Eigen/Dense>
 #include <Eigen/src/Core/util/Constants.h>
 
 #include "common/eigen_tools.hh"
 
 #ifndef T4_MAP_PROXY_H
 #define T4_MAP_PROXY_H
 
 namespace muSpectre {
 
   /**
    * simple adapter function to create a matrix that can be mapped as a tensor
    */
   template <typename T, Dim_t Dim>
   using T4Mat = Eigen::Matrix<T, Dim*Dim, Dim*Dim>;
   template <typename T, Dim_t Dim, bool ConstMap=false>
   using T4MatMap = std::conditional_t<ConstMap,
                                       const Eigen::Map<T4Mat<T, Dim>>,
                                       Eigen::Map<T4Mat<T, Dim>>>;
 
 
   template <typename T4>
-  inline decltype(auto) get(T4&& t4, Dim_t i, Dim_t j, Dim_t k, Dim_t l) {
+  inline auto get(T4&& t4, Dim_t i, Dim_t j, Dim_t k, Dim_t l)
+    -> decltype(t4.coeffRef(i,j)) {
     constexpr Dim_t Dim{EigenCheck::Tensor4Dim(t4)};
     const auto myColStride{
       (t4.colStride() == 1) ? t4.colStride(): t4.colStride()/Dim};
     const auto myRowStride{
       (t4.rowStride() == 1) ? t4.rowStride(): t4.rowStride()/Dim};
     return t4.coeffRef(i * myRowStride + j * myColStride,
                        k * myRowStride + l * myColStride);
 
   }
 
   /* ---------------------------------------------------------------------- */
   /** Proxy class mapping a fourth-order tensor onto a 2D matrix (in
       order to avoid the use of Eigen::Tensor. This class is, however
       byte-compatible with Tensors (i.e., you can map this onto a
       tensor instead of a matrix)
   **/
   template <typename T, Dim_t Dim, bool MapConst=false, bool Symmetric=false,
             int MapOptions=Eigen::Unaligned,
             typename StrideType=Eigen::Stride<0, 0>>
   class T4Map:
     public Eigen::MapBase<T4Map<T, Dim, MapConst, Symmetric,
                                 MapOptions, StrideType>>
   {
     public:
     typedef Eigen::MapBase<T4Map> Base;
     EIGEN_DENSE_PUBLIC_INTERFACE(T4Map);
 
     using matrix_type = Eigen::Matrix<T, Dim*Dim, Dim*Dim>;
     using PlainObjectType =
       std::conditional_t<MapConst,
                          const matrix_type, matrix_type>;
     using ConstType = T4Map<T, Dim, true, Symmetric, MapOptions, StrideType>;
     using Base::colStride;
     using Base::rowStride;
     using Base::IsRowMajor;
     typedef typename Base::PointerType PointerType;
     typedef PointerType PointerArgType;
     using trueScalar = std::conditional_t<MapConst, const Scalar, Scalar>;
     EIGEN_DEVICE_FUNC
     inline PointerType cast_to_pointer_type(PointerArgType ptr) { return ptr; }
 
     EIGEN_DEVICE_FUNC
     inline Eigen::Index innerStride() const
     {
       return StrideType::InnerStrideAtCompileTime != 0 ? m_stride.inner() : 1;
     }
 
     template <class Derived>
     inline T4Map & operator=(const Eigen::MatrixBase<Derived> & other) {
       this->map = other;
       return *this;
     }
 
     EIGEN_DEVICE_FUNC
     inline Eigen::Index outerStride() const
     {
       return StrideType::OuterStrideAtCompileTime != 0 ? m_stride.outer()
         : IsVectorAtCompileTime ? this->size()
         : int(Flags)&Eigen::RowMajorBit ? this->cols()
         : this->rows();
     }
 
     /** Constructor in the fixed-size case.
      *
      * \param dataPtr pointer to the array to map
      * \param stride optional Stride object, passing the strides.
      */
     EIGEN_DEVICE_FUNC
     explicit inline T4Map(PointerArgType dataPtr, const StrideType& stride = StrideType())
       : Base(cast_to_pointer_type(dataPtr)), m_stride(stride),
         map(cast_to_pointer_type(dataPtr))
     {
       PlainObjectType::Base::_check_template_params();
     }
 
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(T4Map);
 
     /** My accessor to mimick tensorial access
      **/
     inline const Scalar& operator()(Dim_t i, Dim_t j, Dim_t k, Dim_t l ) const {
       const auto myColStride{
         (colStride() == 1) ? colStride(): colStride()/Dim};
       const auto myRowStride{
         (rowStride() == 1) ? rowStride(): rowStride()/Dim};
       return this->map.coeff(i * myRowStride + j * myColStride,
                              k * myRowStride + l * myColStride);
     }
 
 
     inline trueScalar& operator()(Dim_t i, Dim_t j, Dim_t k, Dim_t l ) {
       const auto myColStride{
         (colStride() == 1) ? colStride(): colStride()/Dim};
       const auto myRowStride{
         (rowStride() == 1) ? rowStride(): rowStride()/Dim};
       return this->map.coeffRef(i * myRowStride + j * myColStride,
                                 k * myRowStride + l * myColStride);
     }
 
   protected:
     StrideType m_stride;
     Eigen::Map<PlainObjectType> map;
   };
 
 
 }  // muSpectre
 
 namespace Eigen {
   //! forward declarations
   template<typename T> struct traits;
 
   /* ---------------------------------------------------------------------- */
   namespace internal {
     template<typename T, muSpectre::Dim_t Dim, bool MapConst, bool Symmetric,
              int MapOptions, typename StrideType>
     struct traits<muSpectre::T4Map<T, Dim, MapConst, Symmetric,
                                    MapOptions, StrideType> >
       : public traits<Matrix<T, Dim*Dim, Dim*Dim>>
     {
       using PlainObjectType = Matrix<T, Dim*Dim, Dim*Dim>;
       typedef traits<PlainObjectType> TraitsBase;
       enum {
         InnerStrideAtCompileTime = StrideType::InnerStrideAtCompileTime == 0
         ? int(PlainObjectType::InnerStrideAtCompileTime)
         : int(StrideType::InnerStrideAtCompileTime),
         OuterStrideAtCompileTime = StrideType::OuterStrideAtCompileTime == 0
         ? int(PlainObjectType::OuterStrideAtCompileTime)
         : int(StrideType::OuterStrideAtCompileTime),
         Alignment = int(MapOptions)&int(AlignedMask),
         Flags0 = TraitsBase::Flags & (~NestByRefBit),
         Flags = is_lvalue<PlainObjectType>::value ? int(Flags0) : (int(Flags0) & ~LvalueBit)
       };
     private:
       enum { Options }; // Expressions don't have Options
     };
   } // namespace internal
 } // namespace Eigen
 
 #endif /* T4_MAP_PROXY_H */
diff --git a/src/common/tensor_algebra.hh b/src/common/tensor_algebra.hh
index f2d3c83..4404214 100644
--- a/src/common/tensor_algebra.hh
+++ b/src/common/tensor_algebra.hh
@@ -1,278 +1,280 @@
 /**
  * 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
  *
  * @section LICENCE
  *
  * Copyright (C) 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #ifndef TENSOR_ALGEBRA_H
 #define TENSOR_ALGEBRA_H
 
 #include <type_traits>
 
 #include <Eigen/Dense>
 #include <unsupported/Eigen/CXX11/Tensor>
 
 
 #include "common/T4_map_proxy.hh"
 #include "common/common.hh"
 #include "common/eigen_tools.hh"
 
 namespace muSpectre {
 
   namespace Tensors {
 
     template<Dim_t dim>
     using Tens2_t = Eigen::TensorFixedSize<Real, Eigen::Sizes<dim, dim>>;
     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 {
       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 {
 
     template<Dim_t dim>
     using Tens2_t = Eigen::Matrix<Real, dim, dim>;
     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::TensorDim(A)};
       static_assert((dim == EigenCheck::TensorDim(B)),
                      "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::TensorDim(A)};
       static_assert((dim == EigenCheck::TensorDim(B)),
                      "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::TensorDim(A)};
       static_assert((dim == EigenCheck::TensorDim(B)),
                      "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::Tensor4Dim(A)};
       static_assert((dim == EigenCheck::TensorDim(B)),
                     "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));
     }
 
 
   }  // Matrices
 }  // muSpectre
 
 
 #endif /* TENSOR_ALGEBRA_H */
diff --git a/src/common/test_goodies.hh b/src/common/test_goodies.hh
index 47e5e68..fbb0e15 100644
--- a/src/common/test_goodies.hh
+++ b/src/common/test_goodies.hh
@@ -1,70 +1,79 @@
 /**
  * file   test_goodies.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   27 Sep 2017
  *
  * @brief  helpers for testing
  *
  * @section LICENCE
  *
  * Copyright (C) 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include <random>
 #include <type_traits>
 
 #ifndef TEST_GOODIES_H
 #define TEST_GOODIES_H
 
 namespace muSpectre {
 
   namespace testGoodies {
 
+    template <Dim_t Dim>
+    struct dimFixture{
+      constexpr static Dim_t dim{Dim};
+    };
+
+    using dimlist = boost::mpl::list<dimFixture<oneD>,
+                                   dimFixture<twoD>,
+                                   dimFixture<threeD>>;
+
     template<typename T>
     class RandRange {
     public:
       RandRange(): rd(), gen(rd()) {}
 
       template <typename dummyT = T>
       std::enable_if_t<std::is_floating_point<dummyT>::value, dummyT>
       randval(T&& lower, T&& upper) {
         static_assert(std::is_same<T, dummyT>::value);
         auto distro = std::uniform_real_distribution<T>(lower, upper);
         return distro(this->gen);
       }
 
       template <typename dummyT = T>
       std::enable_if_t<std::is_integral<dummyT>::value, dummyT>
       randval(T&& lower, T&& upper) {
         static_assert(std::is_same<T, dummyT>::value);
         auto distro = std::uniform_int_distribution<T>(lower, upper);
         return distro(this->gen);
       }
 
     private:
       std::random_device rd;
       std::default_random_engine gen;
     };
 
   }  // testGoodies
 
 }  // muSpectre
 
 #endif /* TEST_GOODIES_H */
diff --git a/tests/test_materials_toolbox.cc b/tests/test_materials_toolbox.cc
index 3a42ad1..49e2b8d 100644
--- a/tests/test_materials_toolbox.cc
+++ b/tests/test_materials_toolbox.cc
@@ -1,172 +1,166 @@
 /**
  * file   test_materials_toolbox.cc
  *
  * @author Till Junge <till.junge@altermail.ch>
  *
  * @date   05 Nov 2017
  *
  * @brief  Tests for the materials toolbox
  *
  * @section LICENCE
  *
  * Copyright (C) 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 #include <boost/mpl/list.hpp>
 
 #include <Eigen/Dense>
 
 #include "tests.hh"
 #include "materials/materials_toolbox.hh"
 #include "common/T4_map_proxy.hh"
 #include "common/tensor_algebra.hh"
 #include "common/test_goodies.hh"
 
 
 namespace muSpectre {
 
   BOOST_AUTO_TEST_SUITE(materials_toolbox)
 
-  template <Dim_t Dim>
-  struct dimFixture{
-    constexpr static Dim_t dim{Dim};
-  };
-
-  using dimlist = boost::mpl::list<dimFixture<oneD>,
-                                   dimFixture<twoD>,
-                                   dimFixture<threeD>>;
-
-  BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_strain_conversion, Fix, dimlist, Fix){
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_strain_conversion, Fix,
+                                   testGoodies::dimlist, Fix){
     constexpr Dim_t dim{Fix::dim};
     using T2 = Eigen::Matrix<Real, dim, dim>;
 
     T2 F;
     F.setRandom();
     auto Eref = .5*(F.transpose()*F-T2::Identity());
 
     auto E_tb = MatTB::convert_strain<StrainMeasure::Gradient,
                                       StrainMeasure::GreenLagrange>
       (Eigen::Map<Eigen::Matrix<Real, dim, dim>>(F.data()));
 
     auto error = (Eref-E_tb).norm();
     BOOST_CHECK_LT(error, tol);
 
     auto F_tb = MatTB::convert_strain<StrainMeasure::Gradient, StrainMeasure::Gradient>(F);
 
     error = (F-F_tb).norm();
     BOOST_CHECK_LT(error, tol);
   }
 
-  BOOST_FIXTURE_TEST_CASE_TEMPLATE(dumb_tensor_mult_test, Fix, dimlist, Fix) {
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(dumb_tensor_mult_test, Fix,
+                                   testGoodies::dimlist, Fix) {
     constexpr Dim_t dim{Fix::dim};
     using T4 = T4Mat<Real, dim>;
     T4 A,B, R1, R2;
     A.setRandom();
     B.setRandom();
     R1 = A*B;
     R2.setZero();
     for (Dim_t i = 0; i < dim; ++i) {
       for (Dim_t j = 0; j < dim; ++j) {
         for (Dim_t a = 0; a < dim; ++a) {
           for (Dim_t b = 0; b < dim; ++b) {
             for (Dim_t k = 0; k < dim; ++k) {
               for (Dim_t l = 0; l < dim; ++l) {
                 get(R2,i,j,k,l) += get(A, i,j,a,b)*get(B, a,b, k,l);
               }
             }
           }
         }
       }
     }
     auto error{(R1-R2).norm()};
     BOOST_CHECK_LT(error, tol);
   }
 
-  BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_PK1_stress, Fix, dimlist, Fix) {
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_PK1_stress, Fix,
+                                   testGoodies::dimlist, Fix) {
     using namespace Matrices;
     constexpr Dim_t dim{Fix::dim};
     using T2 = Eigen::Matrix<Real, dim, dim>;
     using T4 = T4Mat<Real, dim>;
     testGoodies::RandRange<Real> rng;
     T2 F=T2::Identity()*2 ;
     //F.setRandom();
     T2 E_tb = MatTB::convert_strain<StrainMeasure::Gradient, StrainMeasure::GreenLagrange>
       (Eigen::Map<Eigen::Matrix<Real, dim, dim>>(F.data()));
     Real lambda = 3;//rng.randval(1, 2);
     Real mu = 4;//rng.randval(1,2);
     T4 J = Itrac<dim>();
     T2 I = I2<dim>();
     T4 I4 = Isymm<dim>();
     T4 C = lambda*J + 2*mu*I4;
     T2 S = tensmult(C, E_tb);
     T2 Sref = lambda*E_tb.trace()*I + 2*mu*E_tb;
 
     auto error{(Sref-S).norm()};
     BOOST_CHECK_LT(error, tol);
 
     T4 K = outer_under(I,S) + outer_under(F,I)*C*outer_under(F.transpose(),I);
 
     // See Curnier, 2000, "Méthodes numériques en mécanique des solides", p 252
     T4 Kref;
     Real Fkrkr = (F.array()*F.array()).sum();
     T2 Fkmkn = F.transpose()*F;
     T2 Fisjs = F*F.transpose();
     Kref.setZero();
     for (Dim_t i = 0; i < dim; ++i) {
       for (Dim_t j = 0; j < dim; ++j) {
         for (Dim_t m = 0; m < dim; ++m) {
           for (Dim_t n = 0; n < dim; ++n) {
             get(Kref, i, m, j, n) =
               (lambda*((Fkrkr-dim)/2 * I(i,j)*I(m,n) + F(i,m)*F(j,n)) +
                mu * (I(i,j)*Fkmkn(m,n) + Fisjs(i,j)*I(m,n) -
                      I(i,j) *I(m,n) + F(i,n)*F(j,m)));
           }
         }
       }
     }
     error = (Kref-K).norm();
     BOOST_CHECK_LT(error, tol);
 
     T2 P = MatTB::PK1_stress<StressMeasure::PK2, StrainMeasure::GreenLagrange>(F, S);
     T2 Pref = F*S;
     error = (P-Pref).norm();
     BOOST_CHECK_LT(error, tol);
 
     auto && stress_tgt = MatTB::PK1_stress<StressMeasure::PK2, StrainMeasure::GreenLagrange>(F, S, C);
     T2 P_t = std::move(std::get<0>(stress_tgt));
     T4 K_t = std::move(std::get<1>(stress_tgt));
     error = (P_t-Pref).norm();
     BOOST_CHECK_LT(error, tol);
 
     error = (K_t-Kref).norm();
     BOOST_CHECK_LT(error, tol);
 
     auto && stress_tgt_trivial =
       MatTB::PK1_stress<StressMeasure::PK1, StrainMeasure::Gradient>(F, P, K);
     T2 P_u = std::move(std::get<0>(stress_tgt_trivial));
     T4 K_u = std::move(std::get<1>(stress_tgt_trivial));
 
     error = (P_u-Pref).norm();
     BOOST_CHECK_LT(error, tol);
 
     error = (K_u-Kref).norm();
     BOOST_CHECK_LT(error, tol);
   }
 
   BOOST_AUTO_TEST_SUITE_END();
 
 }  // muSpectre
diff --git a/tests/test_tensor_algebra.cc b/tests/test_tensor_algebra.cc
index 18f4cac..550ee61 100644
--- a/tests/test_tensor_algebra.cc
+++ b/tests/test_tensor_algebra.cc
@@ -1,190 +1,304 @@
 /**
  * file   test_tensor_algebra.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   05 Nov 2017
  *
  * @brief  Tests for the tensor algebra functions
  *
  * @section LICENCE
  *
  * Copyright (C) 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include <iomanip>
 
 #include <unsupported/Eigen/CXX11/Tensor>
 
 #include "common/tensor_algebra.hh"
 #include "tests.hh"
+#include "common/test_goodies.hh"
 
 
 namespace muSpectre {
 
   BOOST_AUTO_TEST_SUITE(tensor_algebra)
   auto TerrNorm = [](auto && t){
     return Eigen::Tensor<Real, 0>(t.abs().sum())();
   };
 
   /* ---------------------------------------------------------------------- */
-  BOOST_AUTO_TEST_CASE(outer_product_test) {
+  BOOST_AUTO_TEST_CASE(tensor_outer_product_test) {
     constexpr Dim_t dim{2};
     Eigen::TensorFixedSize<Real, Eigen::Sizes<dim,dim>> A, B;
     // use prime numbers so that every multiple is uniquely identifiable
     A.setValues({{1, 2}, {3 ,7}});
     B.setValues({{11, 13}, {17 ,19}});
 
     Eigen::TensorFixedSize<Real, Eigen::Sizes<dim, dim, dim, dim>> Res1, Res2, Res3;
 
     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) {
             Res1(i, j, k, l) = A(i, j)*B(k, l);
             Res2(i, j, k, l) = A(i, k)*B(j, l);
             Res3(i, j, k, l) = A(i, l)*B(j, k);
           }
         }
       }
     }
 
 
     Real error = TerrNorm(Res1 - Tensors::outer<dim>(A,B));
     BOOST_CHECK_LT(error, tol);
 
     error = TerrNorm(Res2 - Tensors::outer_under<dim>(A, B));
     BOOST_CHECK_LT(error, tol);
 
     error = TerrNorm(Res3 - Tensors::outer_over<dim>(A, B));
     if (error > tol) {
       std::cout << "reference:" << std::endl
                 << Res3 << std::endl;
       std::cout << "result:" << std::endl
                 << Tensors::outer_over<dim>(A, B) << std::endl;
       std::cout << "A:" << std::endl
                 << A << std::endl;
       std::cout << "B" << std::endl
                 << B << std::endl;
       decltype(Res3) tmp = Tensors::outer_over<dim>(A, B);
       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) {
               std::cout << "for (" << i << ", " << j << ", " << k << ", " << l
                         << "), ref: " << std::setw(3)<< Res3(i,j,k,l)
                         << ", res: " <<  std::setw(3)<< tmp(i,j,k,l)
                         << std::endl;
             }
           }
         }
       }
     }
     BOOST_CHECK_LT(error, tol);
 
     error = TerrNorm(Res3 - Tensors::outer<dim>(A, B));
     BOOST_CHECK_GT(error, tol);
   };
 
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(outer_products, Fix,
+                                   testGoodies::dimlist, Fix) {
+    constexpr auto dim{Fix::dim};
+    using T2 = Tensors::Tens2_t<dim>;
+    using M2 = Matrices::Tens2_t<dim>;
+    using Map2 = Eigen::Map<const M2>;
+    using T4 = Tensors::Tens4_t<dim>;
+    using M4 = Matrices::Tens4_t<dim>;
+    using Map4 = Eigen::Map<const M4>;
+    T2 A, B;
+    T4 RT;
+    A.setRandom();
+    B.setRandom();
+    Map2 Amap(A.data());
+    Map2 Bmap(B.data());
+    M2 C, D;
+    M4 RM;
+    C = Amap;
+    D = Bmap;
+
+    auto error =[](const T4& A, const M4& B) {return (B - Map4(A.data())).norm();};
+
+    // Check outer product
+    RT = Tensors::outer<dim>(A, B);
+    RM = Matrices::outer(C, D);
+    BOOST_CHECK_LT(error(RT, RM), tol);
+
+    // Check outer_under product
+    RT = Tensors::outer_under<dim>(A, B);
+    RM = Matrices::outer_under(C, D);
+    BOOST_CHECK_LT(error(RT, RM), tol);
+
+    // Check outer_over product
+    RT = Tensors::outer_over<dim>(A, B);
+    RM = Matrices::outer_over(C, D);
+    BOOST_CHECK_LT(error(RT, RM), tol);
+
+  }
+
   BOOST_AUTO_TEST_CASE(tensor_multiplication) {
     constexpr Dim_t dim{2};
     using Strain_t = Eigen::TensorFixedSize<Real, Eigen::Sizes<dim, dim>>;
     Strain_t A, B;
     A.setValues({{1, 2}, {3 ,7}});
     B.setValues({{11, 13}, {17 ,19}});
     Strain_t FF1 = A*B; // element-wise multiplication
     std::array<Eigen::IndexPair<int>, 1> prod_dims
     {   Eigen::IndexPair<int>{1, 0}};
 
     Strain_t FF2 = A.contract(B, prod_dims); // correct option 1
     Strain_t FF3;
     using Mat_t = Eigen::Map<Eigen::Matrix<Real, dim, dim>>;
     // following only works for evaluated tensors (which already have data())
     Mat_t(FF3.data()) = Mat_t(A.data()) * Mat_t(B.data());
     Strain_t ref;
     ref.setZero();
     for (Dim_t i = 0; i < dim; ++i) {
       for (Dim_t j = 0; j < dim; ++j) {
         for (Dim_t a = 0; a < dim; ++a) {
           ref(i,j) += A(i, a) * B(a, j);
         }
       }
     }
 
 
     using T = Eigen::TensorFixedSize<Real, Eigen::Sizes<2, 2> >;
     constexpr Dim_t my_size =
       T::NumIndices;
     std::cout << "size = " << my_size << std::endl;
 
     using Strain_tw = Eigen::TensorFixedSize<Real, Eigen::Sizes<dim+1, dim+1>>;
     Strain_tw C;
     C.setConstant(100);
     auto a = T::Dimensions::total_size;
     auto b = Strain_tw::Dimensions::total_size;
     std::cout << "a, b = " << a << ", " << b << std::endl;
     //static_assert(!std::is_convertible<Strain_t, Strain_tw>::value,
     //              "Tensors not size-protected");
     if (std::is_convertible<Strain_t, Strain_tw>::value) {
       std::cout << "this is not good, should I abandon Tensors?";
     }
     // this test seems useless. I use to detect if Eigen changed the
     // default tensor product
     Real error = TerrNorm(FF1-ref);
     if (error < tol) {
       std::cout << "A =" << std::endl
                 << A << std::endl;
       std::cout << "B =" << std::endl
                 << B << std::endl;
       std::cout << "FF1 =" << std::endl
                 << FF1 << std::endl;
       std::cout << "ref =" << std::endl
                 << ref << std::endl;
     }
     BOOST_CHECK_GT(error, tol);
 
     error = TerrNorm(FF2-ref);
     if (error > tol) {
       std::cout << "A =" << std::endl
                 << A << std::endl;
       std::cout << "B =" << std::endl
                 << B << std::endl;
       std::cout << "FF2 =" << std::endl
                 << FF2 << std::endl;
       std::cout << "ref =" << std::endl
                 << ref << std::endl;
     }
     BOOST_CHECK_LT(error, tol);
 
     error = TerrNorm(FF3-ref);
     if (error > tol) {
       std::cout << "A =" << std::endl
                 << A << std::endl;
       std::cout << "B =" << std::endl
                 << B << std::endl;
       std::cout << "FF3 =" << std::endl
                 << FF3 << std::endl;
       std::cout << "ref =" << std::endl
                 << ref << std::endl;
     }
     BOOST_CHECK_LT(error, tol);
 
   }
+
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensmult, Fix,
+                                   testGoodies::dimlist, Fix) {
+    using namespace Matrices;
+    constexpr Dim_t dim{Fix::dim};
+    using T4 = Tens4_t<dim>;
+    using T2 = Tens2_t<dim>;
+    using V2 = Eigen::Matrix<Real, dim*dim, 1>;
+    T4 C;
+    C.setRandom();
+    T2 E;
+    E.setRandom();
+    Eigen::Map<const V2> Ev(E.data());
+    T2 R = tensmult(C, E);
+
+    auto error = (Eigen::Map<const V2>(R.data())-C*Ev).norm();
+    BOOST_CHECK_LT(error, tol);
+  }
+
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tracer, Fix,
+                                   testGoodies::dimlist, Fix) {
+    using namespace Matrices;
+    constexpr Dim_t dim{Fix::dim};
+    using T2 = Tens2_t<dim>;
+    auto tracer = Itrac<dim>();
+    T2 F;
+    F.setRandom();
+    auto Ftrac = tensmult(tracer, F);
+    auto error = (Ftrac-F.trace()*F.Identity()).norm();
+    BOOST_CHECK_LT(error, tol);
+  }
+
+
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_identity, Fix,
+                                   testGoodies::dimlist, Fix) {
+    using namespace Matrices;
+    constexpr Dim_t dim{Fix::dim};
+    using T2 = Tens2_t<dim>;
+    auto ident = Iiden<dim>();
+    T2 F;
+    F.setRandom();
+    auto Fiden = tensmult(ident, F);
+    auto error = (Fiden-F).norm();
+    BOOST_CHECK_LT(error, tol);
+  }
+
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_transposer, Fix,
+                                   testGoodies::dimlist, Fix) {
+    using namespace Matrices;
+    constexpr Dim_t dim{Fix::dim};
+    using T2 = Tens2_t<dim>;
+    auto trnst = Itrns<dim>();
+    T2 F;
+    F.setRandom();
+    auto Ftrns = tensmult(trnst, F);
+    auto error = (Ftrns-F.transpose()).norm();
+    BOOST_CHECK_LT(error, tol);
+  }
+
+  BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_symmetriser, Fix,
+                                   testGoodies::dimlist, Fix) {
+    using namespace Matrices;
+    constexpr Dim_t dim{Fix::dim};
+    using T2 = Tens2_t<dim>;
+    auto symmt = Isymm<dim>();
+    T2 F;
+    F.setRandom();
+    auto Fsymm = tensmult(symmt, F);
+    auto error = (Fsymm-.5*(F+F.transpose())).norm();
+    BOOST_CHECK_LT(error, tol);
+    std::cout << " F =" << std::endl << F << std::endl;
+    std::cout << "IF =" << std::endl << Fsymm << std::endl;
+  }
+
   BOOST_AUTO_TEST_SUITE_END();
 
 }  // muSpectre
diff --git a/tests/tests.hh b/tests/tests.hh
index 03aba7d..cddf95e 100644
--- a/tests/tests.hh
+++ b/tests/tests.hh
@@ -1,44 +1,45 @@
 /**
  * file   tests.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   10 May 2017
  *
  * @brief  common defs for tests
  *
  * @section LICENCE
  *
  * Copyright (C) 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "common/common.hh"
 #include <boost/test/unit_test.hpp>
 #include <boost/mpl/list.hpp>
 
 #ifndef TESTS_H
 #define TESTS_H
 
 namespace muSpectre {
 
+
   const Real tol = 1e-14*100; //it's in percent
 
 }  // muSpectre
 
 
 #endif /* TESTS_H */