diff --git a/test/test_common/test_array.cc b/test/test_common/test_array.cc
index d7d0614f9..10baffdca 100644
--- a/test/test_common/test_array.cc
+++ b/test/test_common/test_array.cc
@@ -1,289 +1,289 @@
 /**
  * @file   test_array.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Nov 09 2017
  * @date last modification: Fri Jan 26 2018
  *
  * @brief  Test the arry class
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "test_gtest_utils.hh"
 /* -------------------------------------------------------------------------- */
 #include <aka_array.hh>
 #include <aka_types.hh>
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <memory>
 #include <typeindex>
 #include <typeinfo>
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 namespace {
 
 class NonTrivial {
 public:
   NonTrivial() = default;
   NonTrivial(int a) : a(a){};
 
   bool operator==(const NonTrivial & rhs) { return a == rhs.a; }
   int a{0};
 };
 
 bool operator==(const int & a, const NonTrivial & rhs) { return a == rhs.a; }
 
 std::ostream & operator<<(std::ostream & stream, const NonTrivial & _this) {
   stream << _this.a;
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 using TestTypes = ::testing::Types<Real, UInt, NonTrivial>;
 /* -------------------------------------------------------------------------- */
 
 ::testing::AssertionResult AssertType(const char * /*a_expr*/,
                                       const char * /*b_expr*/,
                                       const std::type_info & a,
                                       const std::type_info & b) {
   if (std::type_index(a) == std::type_index(b))
     return ::testing::AssertionSuccess();
 
   return ::testing::AssertionFailure()
          << debug::demangle(a.name()) << " != " << debug::demangle(b.name())
          << ") are different";
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T> class ArrayConstructor : public ::testing::Test {
 protected:
   using type = T;
 
   void SetUp() override { type_str = debug::demangle(typeid(T).name()); }
 
   template <typename... P> decltype(auto) construct(P &&... params) {
     return std::make_unique<Array<T>>(std::forward<P>(params)...);
   }
 
 protected:
   std::string type_str;
 };
 
-TYPED_TEST_SUITE(ArrayConstructor, TestTypes);
+TYPED_TEST_SUITE(ArrayConstructor, TestTypes, );
 
 TYPED_TEST(ArrayConstructor, ConstructDefault1) {
   auto array = this->construct();
   EXPECT_EQ(0, array->size());
   EXPECT_EQ(1, array->getNbComponent());
   EXPECT_STREQ("", array->getID().c_str());
 }
 
 TYPED_TEST(ArrayConstructor, ConstructDefault2) {
   auto array = this->construct(1000);
   EXPECT_EQ(1000, array->size());
   EXPECT_EQ(1, array->getNbComponent());
   EXPECT_STREQ("", array->getID().c_str());
 }
 
 TYPED_TEST(ArrayConstructor, ConstructDefault3) {
   auto array = this->construct(1000, 10);
   EXPECT_EQ(1000, array->size());
   EXPECT_EQ(10, array->getNbComponent());
   EXPECT_STREQ("", array->getID().c_str());
 }
 
 TYPED_TEST(ArrayConstructor, ConstructDefault4) {
   auto array = this->construct(1000, 10, "test");
   EXPECT_EQ(1000, array->size());
   EXPECT_EQ(10, array->getNbComponent());
   EXPECT_STREQ("test", array->getID().c_str());
 }
 
 TYPED_TEST(ArrayConstructor, ConstructDefault5) {
   auto array = this->construct(1000, 10, 1);
   EXPECT_EQ(1000, array->size());
   EXPECT_EQ(10, array->getNbComponent());
   EXPECT_EQ(1, array->operator()(10, 6));
   EXPECT_STREQ("", array->getID().c_str());
 }
 
 // TYPED_TEST(ArrayConstructor, ConstructDefault6) {
 //   typename TestFixture::type defaultv[2] = {0, 1};
 
 //   auto array = this->construct(1000, 2, defaultv);
 //   EXPECT_EQ(1000, array->size());
 //   EXPECT_EQ(2, array->getNbComponent());
 //   EXPECT_EQ(1, array->operator()(10, 1));
 //   EXPECT_EQ(0, array->operator()(603, 0));
 //   EXPECT_STREQ("", array->getID().c_str());
 // }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> class ArrayFixture : public ArrayConstructor<T> {
 public:
   void SetUp() override {
     ArrayConstructor<T>::SetUp();
     array = this->construct(1000, 10);
   }
 
   void TearDown() override { array.reset(nullptr); }
 
 protected:
   std::unique_ptr<Array<T>> array;
 };
 
-TYPED_TEST_SUITE(ArrayFixture, TestTypes);
+TYPED_TEST_SUITE(ArrayFixture, TestTypes, );
 
 TYPED_TEST(ArrayFixture, Copy) {
   Array<typename TestFixture::type> copy(*this->array);
 
   EXPECT_EQ(1000, copy.size());
   EXPECT_EQ(10, copy.getNbComponent());
   EXPECT_NE(this->array->storage(), copy.storage());
 }
 
 TYPED_TEST(ArrayFixture, Set) {
   auto & arr = *(this->array);
   arr.set(12);
   EXPECT_EQ(12, arr(156, 5));
   EXPECT_EQ(12, arr(520, 7));
   EXPECT_EQ(12, arr(999, 9));
 }
 
 TYPED_TEST(ArrayFixture, Resize) {
   auto & arr = *(this->array);
 
   auto * ptr = arr.storage();
 
   arr.resize(0);
   EXPECT_EQ(0, arr.size());
   EXPECT_TRUE(arr.storage() == nullptr or arr.storage() == ptr);
   EXPECT_LE(0, arr.getAllocatedSize());
 
   arr.resize(3000);
   EXPECT_EQ(3000, arr.size());
   EXPECT_LE(3000, arr.getAllocatedSize());
 
   ptr = arr.storage();
 
   arr.resize(0);
   EXPECT_EQ(0, arr.size());
   EXPECT_TRUE(arr.storage() == nullptr or arr.storage() == ptr);
   EXPECT_LE(0, arr.getAllocatedSize());
 }
 
 TYPED_TEST(ArrayFixture, PushBack) {
   auto & arr = *(this->array);
 
   auto * ptr = arr.storage();
 
   arr.resize(0);
   EXPECT_EQ(0, arr.size());
   EXPECT_TRUE(arr.storage() == nullptr or arr.storage() == ptr);
   EXPECT_LE(0, arr.getAllocatedSize());
 
   arr.resize(3000);
   EXPECT_EQ(3000, arr.size());
   EXPECT_LE(3000, arr.getAllocatedSize());
 
   ptr = arr.storage();
 
   arr.resize(0);
   EXPECT_EQ(0, arr.size());
   EXPECT_TRUE(arr.storage() == nullptr or arr.storage() == ptr);
   EXPECT_LE(0, arr.getAllocatedSize());
 }
 
 TYPED_TEST(ArrayFixture, ViewVector) {
   auto && view = make_view(*this->array, 10);
   EXPECT_NO_THROW(view.begin());
   {
     auto it = view.begin();
     EXPECT_EQ(10, it->size());
     EXPECT_PRED_FORMAT2(AssertType, typeid(*it),
                         typeid(Vector<typename TestFixture::type>));
     EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]),
                         typeid(VectorProxy<typename TestFixture::type>));
   }
 }
 
 TYPED_TEST(ArrayFixture, ViewMatrix) {
   {
     auto && view = make_view(*this->array, 2, 5);
 
     EXPECT_NO_THROW(view.begin());
     {
       auto it = view.begin();
       EXPECT_EQ(10, it->size());
       EXPECT_EQ(2, it->size(0));
       EXPECT_EQ(5, it->size(1));
 
       EXPECT_PRED_FORMAT2(AssertType, typeid(*it),
                           typeid(Matrix<typename TestFixture::type>));
       EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]),
                           typeid(MatrixProxy<typename TestFixture::type>));
     }
   }
 }
 
 TYPED_TEST(ArrayFixture, ViewVectorWrong) {
   auto && view = make_view(*this->array, 11);
   EXPECT_THROW(view.begin(), debug::ArrayException);
 }
 
 TYPED_TEST(ArrayFixture, ViewMatrixWrong) {
   auto && view = make_view(*this->array, 3, 7);
   EXPECT_THROW(view.begin(), debug::ArrayException);
 }
 
 TYPED_TEST(ArrayFixture, ViewMatrixIter) {
   std::size_t count = 0;
   for (auto && mat : make_view(*this->array, 10, 10)) {
     EXPECT_EQ(100, mat.size());
     EXPECT_EQ(10, mat.size(0));
     EXPECT_EQ(10, mat.size(1));
     EXPECT_PRED_FORMAT2(AssertType, typeid(mat),
                         typeid(Matrix<typename TestFixture::type>));
 
     ++count;
   }
 
   EXPECT_EQ(100, count);
 }
 
 TYPED_TEST(ArrayFixture, ConstViewVector) {
   const auto & carray = *this->array;
   auto && view = make_view(carray, 10);
   EXPECT_NO_THROW(view.begin());
   {
     auto it = view.begin();
     EXPECT_EQ(10, it->size());
     EXPECT_PRED_FORMAT2(AssertType, typeid(*it),
                         typeid(Vector<typename TestFixture::type>));
     EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]),
                         typeid(VectorProxy<typename TestFixture::type>));
   }
 }
 
 } // namespace
diff --git a/test/test_common/test_voigt_helper.cc b/test/test_common/test_voigt_helper.cc
index 34f4a33b6..33a9440b9 100644
--- a/test/test_common/test_voigt_helper.cc
+++ b/test/test_common/test_voigt_helper.cc
@@ -1,160 +1,160 @@
 /**
  * @file   test_voigt_helper.cc
  *
  * @author Nicolas Richart
  *
  * @date creation  mer nov 13 2019
  *
  * @brief unit tests for VoigtHelper
  *
  *
  * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "test_gtest_utils.hh"
 /* -------------------------------------------------------------------------- */
 #include <aka_voigthelper.hh>
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <unordered_map>
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 template <class Dim_v> class VoigtHelperFixture : public ::testing::Test {
 protected:
   using voigt_h = VoigtHelper<Dim_v::value>;
   constexpr static UInt dim = Dim_v::value;
 
   VoigtHelperFixture() {
     switch (this->dim) {
     case 1: {
       indices.push_back({0, 0});
       matrix = Matrix<Real>{{10}};
       vector = Vector<Real>{10};
       vector_factor = Vector<Real>{10};
       break;
     }
     case 2: {
       indices.push_back({0, 0});
       indices.push_back({1, 1});
       indices.push_back({0, 1});
 
       matrix = Matrix<Real>{{10, 33}, {0, 56}};
       vector = Vector<Real>{10, 56, 33};
       vector_factor = Vector<Real>{10, 56, 2 * 33};
       break;
     }
     case 3: {
       indices.push_back({0, 0});
       indices.push_back({1, 1});
       indices.push_back({2, 2});
       indices.push_back({1, 2});
       indices.push_back({0, 2});
       indices.push_back({0, 1});
 
       matrix = Matrix<Real>{{10, 33, 20}, {0, 56, 27}, {0, 0, 98}};
       vector = Vector<Real>{10, 56, 98, 27, 20, 33};
       vector_factor = Vector<Real>{10, 56, 98, 2 * 27, 2 * 20, 2 * 33};
       break;
     }
     }
   }
 
   void SetUp() override {}
 
   std::vector<std::pair<UInt, UInt>> indices;
   Matrix<Real> matrix;
   Vector<Real> vector;
   Vector<Real> vector_factor;
 };
 
 template <UInt dim>
 using spatial_dimension_t = std::integral_constant<UInt, dim>;
 
 using TestTypes =
     ::testing::Types<spatial_dimension_t<1>, spatial_dimension_t<2>,
                      spatial_dimension_t<3>>;
-TYPED_TEST_SUITE(VoigtHelperFixture, TestTypes);
+TYPED_TEST_SUITE(VoigtHelperFixture, TestTypes, );
 
 TYPED_TEST(VoigtHelperFixture, Size) {
   using voigt_h = typename TestFixture::voigt_h;
   switch (this->dim) {
   case 1:
     EXPECT_EQ(voigt_h::size, 1);
     break;
   case 2:
     EXPECT_EQ(voigt_h::size, 3);
     break;
   case 3:
     EXPECT_EQ(voigt_h::size, 6);
     break;
   }
 }
 
 TYPED_TEST(VoigtHelperFixture, Indicies) {
   using voigt_h = typename TestFixture::voigt_h;
 
   for (UInt I = 0; I < voigt_h::size; ++I) {
     EXPECT_EQ(this->indices[I].first, voigt_h::vec[I][0]);
     EXPECT_EQ(this->indices[I].second, voigt_h::vec[I][1]);
   }
 }
 
 TYPED_TEST(VoigtHelperFixture, Factors) {
   using voigt_h = typename TestFixture::voigt_h;
   for (UInt I = 0; I < voigt_h::size; ++I) {
     if (I < this->dim) {
       EXPECT_EQ(voigt_h::factors[I], 1);
     } else {
       EXPECT_EQ(voigt_h::factors[I], 2);
     }
   }
 }
 
 TYPED_TEST(VoigtHelperFixture, MatrixToVoight) {
   using voigt_h = typename TestFixture::voigt_h;
 
   auto voigt = voigt_h::matrixToVoigt(this->matrix);
 
   for (UInt I = 0; I < voigt_h::size; ++I) {
     EXPECT_EQ(voigt(I), this->vector(I));
   }
 }
 
 TYPED_TEST(VoigtHelperFixture, MatrixToVoightFactors) {
   using voigt_h = typename TestFixture::voigt_h;
 
   auto voigt = voigt_h::matrixToVoigtWithFactors(this->matrix);
 
   for (UInt I = 0; I < voigt_h::size; ++I) {
     EXPECT_EQ(voigt(I), this->vector_factor(I));
   }
 }
 
 TYPED_TEST(VoigtHelperFixture, VoightToMatrix) {
   using voigt_h = typename TestFixture::voigt_h;
 
   auto matrix = voigt_h::voigtToMatrix(this->vector);
 
   for (UInt i = 0; i < this->dim; ++i) {
     for (UInt j = 0; j < this->dim; ++j) {
       EXPECT_EQ(matrix(i, j), this->matrix(std::min(i, j), std::max(i, j)));
     }
   }
 }
diff --git a/test/test_fe_engine/test_fe_engine_fixture.hh b/test/test_fe_engine/test_fe_engine_fixture.hh
index 700d72b85..c3662b687 100644
--- a/test/test_fe_engine/test_fe_engine_fixture.hh
+++ b/test/test_fe_engine/test_fe_engine_fixture.hh
@@ -1,111 +1,111 @@
 /**
  * @file   test_fe_engine_fixture.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Nov 14 2017
  * @date last modification: Mon Feb 19 2018
  *
  * @brief  Fixture for feengine tests
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "test_gtest_utils.hh"
 /* -------------------------------------------------------------------------- */
 #include <element_class.hh>
 #include <fe_engine.hh>
 #include <integrator_gauss.hh>
 #include <shape_lagrange.hh>
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_TEST_FE_ENGINE_FIXTURE_HH_
 #define AKANTU_TEST_FE_ENGINE_FIXTURE_HH_
 
 using namespace akantu;
 
 /// Generic class for FEEngine tests
 template <typename type_, template <ElementKind> class shape_t,
           ElementKind kind = _ek_regular>
 class TestFEMBaseFixture : public ::testing::Test {
 public:
   static constexpr const ElementType type = type_::value;
   static constexpr const size_t dim = ElementClass<type>::getSpatialDimension();
   using FEM = FEEngineTemplate<IntegratorGauss, shape_t, kind>;
 
   /// Setup reads mesh corresponding to element type and initializes an FEEngine
   void SetUp() override {
     const auto dim = this->dim;
     mesh = std::make_unique<Mesh>(dim);
 
     std::stringstream meshfilename;
     meshfilename << type << ".msh";
     this->readMesh(meshfilename.str());
 
     lower = mesh->getLowerBounds();
     upper = mesh->getUpperBounds();
 
     nb_element = this->mesh->getNbElement(type);
 
     fem = std::make_unique<FEM>(*mesh, dim, "my_fem");
     nb_quadrature_points_total =
         GaussIntegrationElement<type>::getNbQuadraturePoints() * nb_element;
 
     SCOPED_TRACE(std::to_string(type));
   }
 
   void TearDown() override {
     fem.reset(nullptr);
     mesh.reset(nullptr);
   }
 
   /// Should be reimplemented if further treatment of the mesh is needed
   virtual void readMesh(std::string file_name) { mesh->read(file_name); }
 
 protected:
   std::unique_ptr<FEM> fem;
   std::unique_ptr<Mesh> mesh;
   UInt nb_element;
   UInt nb_quadrature_points_total;
   Vector<Real> lower;
   Vector<Real> upper;
 };
 
 template <typename type_, template <ElementKind> class shape_t,
           ElementKind kind>
 constexpr const ElementType TestFEMBaseFixture<type_, shape_t, kind>::type;
 
 template <typename type_, template <ElementKind> class shape_t,
           ElementKind kind>
 constexpr const size_t TestFEMBaseFixture<type_, shape_t, kind>::dim;
 
 /* -------------------------------------------------------------------------- */
 /// Base class for test with Lagrange FEEngine and regular elements
 template <typename type_>
 using TestFEMFixture = TestFEMBaseFixture<type_, ShapeLagrange, _ek_regular>;
 
 /* -------------------------------------------------------------------------- */
 
 using fe_engine_types = gtest_list_t<TestElementTypes>;
 
-TYPED_TEST_SUITE(TestFEMFixture, fe_engine_types);
+TYPED_TEST_SUITE(TestFEMFixture, fe_engine_types, );
 
 #endif /* AKANTU_TEST_FE_ENGINE_FIXTURE_HH_ */
diff --git a/test/test_fe_engine/test_fe_engine_gauss_integration.cc b/test/test_fe_engine/test_fe_engine_gauss_integration.cc
index 8391175f3..e9ea2dde8 100644
--- a/test/test_fe_engine/test_fe_engine_gauss_integration.cc
+++ b/test/test_fe_engine/test_fe_engine_gauss_integration.cc
@@ -1,153 +1,153 @@
 /**
  * @file   test_fe_engine_gauss_integration.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue May 24 2016
  * @date last modification: Mon Feb 19 2018
  *
  * @brief  test integration on elements, this test consider that mesh is a cube
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "test_fe_engine_fixture.hh"
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <iostream>
 /* -------------------------------------------------------------------------- */
 using namespace akantu;
 
 namespace {
 /* -------------------------------------------------------------------------- */
 template <size_t t> using degree_t = std::integral_constant<size_t, t>;
 
 /* -------------------------------------------------------------------------- */
 using TestDegreeTypes = std::tuple<degree_t<0>, degree_t<1>, degree_t<2>,
                                    degree_t<3>, degree_t<4>, degree_t<5>>;
 
 std::array<Polynomial<5>, 3> global_polys{
     {{0.40062394, 0.13703225, 0.51731446, 0.87830084, 0.5410543, 0.71842292},
      {0.41861835, 0.11080576, 0.49874043, 0.49077504, 0.85073835, 0.66259755},
      {0.92620845, 0.7503478, 0.62962232, 0.31662719, 0.64069644, 0.30878135}}};
 
 template <typename T>
 class TestGaussIntegrationFixture
     : public TestFEMFixture<std::tuple_element_t<0, T>> {
 protected:
   using parent = TestFEMFixture<std::tuple_element_t<0, T>>;
   static constexpr size_t degree{std::tuple_element_t<1, T>::value};
 
 public:
   TestGaussIntegrationFixture() : integration_points_pos(0, parent::dim) {}
 
   void SetUp() override {
     parent::SetUp();
     this->fem->initShapeFunctions();
 
     auto integration_points =
              this->fem->getIntegrator().template getIntegrationPoints <
              parent::type,
          degree == 0 ? 1 : degree > ();
 
     nb_integration_points = integration_points.cols();
 
     auto shapes_size = ElementClass<parent::type>::getShapeSize();
     Array<Real> shapes(0, shapes_size);
     this->fem->getShapeFunctions()
         .template computeShapesOnIntegrationPoints<parent::type>(
             this->mesh->getNodes(), integration_points, shapes, _not_ghost);
 
     auto vect_size = this->nb_integration_points * this->nb_element;
     integration_points_pos.resize(vect_size);
     this->fem->getShapeFunctions()
         .template interpolateOnIntegrationPoints<parent::type>(
             this->mesh->getNodes(), integration_points_pos, this->dim, shapes);
 
     for (size_t d = 0; d < this->dim; ++d) {
       polys[d] = global_polys[d].extract(degree);
     }
   }
 
   void testIntegrate() {
     std::stringstream sstr;
     sstr << this->type << ":" << this->degree;
     SCOPED_TRACE(sstr.str().c_str());
 
     auto vect_size = this->nb_integration_points * this->nb_element;
     Array<Real> polynomial(vect_size);
     size_t dim = parent::dim;
 
     for (size_t d = 0; d < dim; ++d) {
       auto poly = this->polys[d];
       for (auto && pair :
            zip(polynomial, make_view(this->integration_points_pos, dim))) {
         auto && p = std::get<0>(pair);
         auto & x = std::get<1>(pair);
         p = poly(x(d));
       }
 
       auto res =
           this->fem->getIntegrator()
               .template integrate<parent::type, (degree == 0 ? 1 : degree)>(
                   polynomial);
       auto expect = poly.integrate(this->lower(d), this->upper(d));
 
       for (size_t o = 0; o < dim; ++o) {
         if (o == d)
           continue;
         expect *= this->upper(d) - this->lower(d);
       }
 
       EXPECT_NEAR(expect, res, 5e-14);
     }
   }
 
 protected:
   UInt nb_integration_points;
   std::array<Array<Real>, parent::dim> polynomial;
   Array<Real> integration_points_pos;
   std::array<Polynomial<5>, 3> polys;
 };
 
 template <typename T> constexpr size_t TestGaussIntegrationFixture<T>::degree;
 
 /* -------------------------------------------------------------------------- */
 /* Tests                                                                      */
 /* -------------------------------------------------------------------------- */
 TYPED_TEST_SUITE_P(TestGaussIntegrationFixture);
 
 TYPED_TEST_P(TestGaussIntegrationFixture, ArbitraryOrder) {
   this->testIntegrate();
 }
 
 REGISTER_TYPED_TEST_SUITE_P(TestGaussIntegrationFixture, ArbitraryOrder);
 
 using TestTypes = gtest_list_t<
     tuple_split_t<50, cross_product_t<TestElementTypes, TestDegreeTypes>>>;
 
-INSTANTIATE_TYPED_TEST_SUITE_P(Split1, TestGaussIntegrationFixture, TestTypes);
+INSTANTIATE_TYPED_TEST_SUITE_P(Split1, TestGaussIntegrationFixture, TestTypes, );
 
 using TestTypesTail = gtest_list_t<
     tuple_split_tail_t<50, cross_product_t<TestElementTypes, TestDegreeTypes>>>;
 
 INSTANTIATE_TYPED_TEST_SUITE_P(Split2, TestGaussIntegrationFixture,
-                               TestTypesTail);
+                               TestTypesTail, );
 } // namespace
diff --git a/test/test_fe_engine/test_fe_engine_precomputation.cc b/test/test_fe_engine/test_fe_engine_precomputation.cc
index 88ba8f632..e7c6a5f35 100644
--- a/test/test_fe_engine/test_fe_engine_precomputation.cc
+++ b/test/test_fe_engine/test_fe_engine_precomputation.cc
@@ -1,113 +1,113 @@
 /**
  * @file   test_fe_engine_precomputation.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Mon Jun 14 2010
  * @date last modification: Mon Feb 19 2018
  *
  * @brief  test of the fem class
  *
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "py_aka_array.hh"
 #include "test_fe_engine_fixture.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/embed.h>
 #include <pybind11/numpy.h>
 /* -------------------------------------------------------------------------- */
 using namespace akantu;
 
 namespace py = pybind11;
 using namespace py::literals;
 
 template <class T> decltype(auto) make_proxy(Array<T> & array) {
   return detail::ArrayProxy<T>(array);
 }
 
 template <typename type_>
 class TestFEMPyFixture : public TestFEMFixture<type_> {
   using parent = TestFEMFixture<type_>;
 
 public:
   void SetUp() override {
     parent::SetUp();
 
     const auto & connectivities = this->mesh->getConnectivity(this->type);
     const auto & nodes = this->mesh->getNodes().begin(this->dim);
     coordinates = std::make_unique<Array<Real>>(
         connectivities.size(), connectivities.getNbComponent() * this->dim);
 
     for (auto && tuple :
          zip(make_view(connectivities, connectivities.getNbComponent()),
              make_view(*coordinates, this->dim,
                        connectivities.getNbComponent()))) {
       const auto & conn = std::get<0>(tuple);
       const auto & X = std::get<1>(tuple);
       for (auto s : arange(conn.size())) {
         Vector<Real>(X(s)) = Vector<Real>(nodes[conn(s)]);
       }
     }
   }
 
   void TearDown() override {
     parent::TearDown();
     coordinates.reset(nullptr);
   }
 
 protected:
   std::unique_ptr<Array<Real>> coordinates;
 };
 
-TYPED_TEST_SUITE(TestFEMPyFixture, fe_engine_types);
+TYPED_TEST_SUITE(TestFEMPyFixture, fe_engine_types, );
 
 TYPED_TEST(TestFEMPyFixture, Precompute) {
   SCOPED_TRACE(std::to_string(this->type));
   this->fem->initShapeFunctions();
   const auto & N = this->fem->getShapeFunctions().getShapes(this->type);
   const auto & B =
       this->fem->getShapeFunctions().getShapesDerivatives(this->type);
   const auto & j = this->fem->getIntegrator().getJacobians(this->type);
 
   // Array<Real> ref_N(this->nb_quadrature_points_total, N.getNbComponent());
   // Array<Real> ref_B(this->nb_quadrature_points_total, B.getNbComponent());
   Array<Real> ref_j(this->nb_quadrature_points_total, j.getNbComponent());
   auto ref_N(N);
   auto ref_B(B);
   py::module py_engine = py::module::import("py_engine");
   auto py_shape = py_engine.attr("Shapes")(py::str(std::to_string(this->type)));
   auto kwargs = py::dict("N"_a = ref_N, "B"_a = ref_B, "j"_a = ref_j,
                          "X"_a = *this->coordinates,
                          "Q"_a = this->fem->getIntegrationPoints(this->type));
 
   auto ret = py_shape.attr("precompute")(**kwargs);
   auto check = [&](auto & ref_A, auto & A, const auto & id) {
     SCOPED_TRACE(std::to_string(this->type) + " " + id);
     for (auto && n : zip(make_view(ref_A, ref_A.getNbComponent()),
                          make_view(A, A.getNbComponent()))) {
       auto diff = (std::get<0>(n) - std::get<1>(n)).template norm<L_inf>();
       EXPECT_NEAR(0., diff, 1e-10);
     }
   };
   check(ref_N, N, "N");
   check(ref_B, B, "B");
   check(ref_j, j, "j");
 }
diff --git a/test/test_model/patch_tests/patch_test_linear_heat_transfer_fixture.hh b/test/test_model/patch_tests/patch_test_linear_heat_transfer_fixture.hh
index 449544034..f1ea709aa 100644
--- a/test/test_model/patch_tests/patch_test_linear_heat_transfer_fixture.hh
+++ b/test/test_model/patch_tests/patch_test_linear_heat_transfer_fixture.hh
@@ -1,76 +1,76 @@
 /**
  * @file   patch_test_linear_heat_transfer_fixture.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Jan 30 2018
  * @date last modification: Wed Jan 31 2018
  *
  * @brief  HeatTransfer patch tests fixture
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "heat_transfer_model.hh"
 /* -------------------------------------------------------------------------- */
 #include "patch_test_linear_fixture.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_PATCH_TEST_LINEAR_HEAT_TRANSFER_FIXTURE_HH_
 #define AKANTU_PATCH_TEST_LINEAR_HEAT_TRANSFER_FIXTURE_HH_
 
 /* -------------------------------------------------------------------------- */
 template <typename type>
 class TestPatchTestHTMLinear
     : public TestPatchTestLinear<type, HeatTransferModel> {
   using parent = TestPatchTestLinear<type, HeatTransferModel>;
 
 public:
   void applyBC() override {
     parent::applyBC();
     auto & temperature = this->model->getTemperature();
     this->applyBConDOFs(temperature);
   }
 
   void initModel(const AnalysisMethod & method,
                  const std::string & material_file) override {
     TestPatchTestLinear<type, HeatTransferModel>::initModel(method,
                                                             material_file);
     if (method != _static)
       this->model->setTimeStep(0.5 * this->model->getStableTimeStep());
   }
 
   void checkAll() {
     auto & temperature = this->model->getTemperature();
     Matrix<Real> C = this->model->get("conductivity");
     this->checkDOFs(temperature);
     this->checkGradient(this->model->getTemperatureGradient(this->type),
                         temperature);
     this->checkResults(
         [&](const Matrix<Real> & grad_T) { return C * grad_T.transpose(); },
         this->model->getKgradT(this->type), temperature);
   }
 };
 
 using htm_types = gtest_list_t<TestElementTypes>;
 
-TYPED_TEST_SUITE(TestPatchTestHTMLinear, htm_types);
+TYPED_TEST_SUITE(TestPatchTestHTMLinear, htm_types, );
 
 #endif /* AKANTU_PATCH_TEST_LINEAR_HEAT_TRANSFER_FIXTURE_HH_ */
diff --git a/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.hh b/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.hh
index 1ac362d9c..2cd9e46aa 100644
--- a/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.hh
+++ b/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.hh
@@ -1,153 +1,153 @@
 /**
  * @file   patch_test_linear_solid_mechanics_fixture.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Jan 30 2018
  *
  * @brief  SolidMechanics patch tests fixture
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "patch_test_linear_fixture.hh"
 #include "solid_mechanics_model.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_PATCH_TEST_LINEAR_SOLID_MECHANICS_FIXTURE_HH_
 #define AKANTU_PATCH_TEST_LINEAR_SOLID_MECHANICS_FIXTURE_HH_
 
 /* -------------------------------------------------------------------------- */
 template <typename tuple_>
 class TestPatchTestSMMLinear
     : public TestPatchTestLinear<std::tuple_element_t<0, tuple_>,
                                  SolidMechanicsModel> {
   using parent =
       TestPatchTestLinear<std::tuple_element_t<0, tuple_>, SolidMechanicsModel>;
 
 public:
   static constexpr bool plane_strain = std::tuple_element_t<1, tuple_>::value;
 
   void applyBC() override {
     parent::applyBC();
     auto & displacement = this->model->getDisplacement();
     this->applyBConDOFs(displacement);
   }
 
   void checkForces() {
     auto & mat = this->model->getMaterial(0);
     auto & internal_forces = this->model->getInternalForce();
     auto & external_forces = this->model->getExternalForce();
     auto dim = this->dim;
 
     Matrix<Real> sigma =
         make_view(mat.getStress(this->type), dim, dim).begin()[0];
 
     external_forces.zero();
     if (dim > 1) {
       for (auto & eg : this->mesh->iterateElementGroups()) {
         this->model->applyBC(BC::Neumann::FromHigherDim(sigma), eg.getName());
       }
     } else {
       external_forces(0) = -sigma(0, 0);
       external_forces(1) = sigma(0, 0);
     }
 
     Real force_norm_inf = -std::numeric_limits<Real>::max();
 
     Vector<Real> total_force(dim);
     total_force.zero();
 
     for (auto && f : make_view(internal_forces, dim)) {
       total_force += f;
       force_norm_inf = std::max(force_norm_inf, f.template norm<L_inf>());
     }
 
     EXPECT_NEAR(0, total_force.template norm<L_inf>() / force_norm_inf, 1e-9);
 
     for (auto && tuple : zip(make_view(internal_forces, dim),
                              make_view(external_forces, dim))) {
       auto && f_int = std::get<0>(tuple);
       auto && f_ext = std::get<1>(tuple);
       auto f = f_int + f_ext;
       EXPECT_NEAR(0, f.template norm<L_inf>() / force_norm_inf, 1e-9);
     }
   }
 
   void checkAll() {
     auto & displacement = this->model->getDisplacement();
     auto & mat = this->model->getMaterial(0);
 
     this->checkDOFs(displacement);
     this->checkGradient(mat.getGradU(this->type), displacement);
     this->checkResults(
         [&](const Matrix<Real> & pstrain) {
           Real nu = this->model->getMaterial(0).get("nu");
           Real E = this->model->getMaterial(0).get("E");
 
           auto strain = (pstrain + pstrain.transpose()) / 2.;
           auto trace = strain.trace();
 
           auto lambda = nu * E / ((1 + nu) * (1 - 2 * nu));
           auto mu = E / (2 * (1 + nu));
 
           if (not this->plane_strain) {
             lambda = nu * E / (1 - nu * nu);
           }
 
           decltype(strain) stress(this->dim, this->dim);
 
           if (this->dim == 1) {
             stress(0, 0) = E * strain(0, 0);
           } else {
             for (UInt i = 0; i < this->dim; ++i)
               for (UInt j = 0; j < this->dim; ++j)
                 stress(i, j) =
                     (i == j) * lambda * trace + 2 * mu * strain(i, j);
           }
 
           return stress;
         },
         mat.getStress(this->type), displacement);
     this->checkForces();
   }
 };
 
 template <typename tuple_>
 constexpr bool TestPatchTestSMMLinear<tuple_>::plane_strain;
 
 template <typename T> struct invalid_plan_stress : std::true_type {};
 template <typename type, typename bool_c>
 struct invalid_plan_stress<std::tuple<type, bool_c>>
     : aka::bool_constant<ElementClass<type::value>::getSpatialDimension() !=
                              2 and
                          not bool_c::value> {};
 
 using true_false =
     std::tuple<aka::bool_constant<true>, aka::bool_constant<false>>;
 
 template <typename T> using valid_types = aka::negation<invalid_plan_stress<T>>;
 
 using model_types = gtest_list_t<
     tuple_filter_t<valid_types, cross_product_t<TestElementTypes, true_false>>>;
 
-TYPED_TEST_SUITE(TestPatchTestSMMLinear, model_types);
+TYPED_TEST_SUITE(TestPatchTestSMMLinear, model_types, );
 
 #endif /* AKANTU_PATCH_TEST_LINEAR_SOLID_MECHANICS_FIXTURE_HH_ */
diff --git a/test/test_model/patch_tests/test_lumped_mass.cc b/test/test_model/patch_tests/test_lumped_mass.cc
index 33099188a..f149fcaa1 100644
--- a/test/test_model/patch_tests/test_lumped_mass.cc
+++ b/test/test_model/patch_tests/test_lumped_mass.cc
@@ -1,101 +1,101 @@
 /**
  * @file   test_lumped_mass.cc
  *
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Dec 05 2017
  * @date last modification: Tue Jan 30 2018
  *
  * @brief  test the lumping of the mass matrix
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "solid_mechanics_model.hh"
 #include "test_gtest_utils.hh"
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <tuple>
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 template <typename tuple_>
 class TestLumpedMassesFixture : public ::testing::Test {
 public:
   static constexpr ElementType type = tuple_::value;
   static constexpr size_t dim = ElementClass<type>::getSpatialDimension();
 
   void SetUp() override {
     debug::setDebugLevel(dblError);
     getStaticParser().parse("material_lumped.dat");
 
     std::stringstream element_type;
     element_type << type;
 
     mesh = std::make_unique<Mesh>(dim);
     mesh->read(element_type.str() + ".msh");
 
     SCOPED_TRACE(element_type.str().c_str());
 
     model = std::make_unique<SolidMechanicsModel>(*mesh);
 
     model->initFull(_analysis_method = _explicit_lumped_mass);
   }
 
   void TearDown() override {
     model.reset(nullptr);
     mesh.reset(nullptr);
   }
 
 protected:
   std::unique_ptr<Mesh> mesh;
   std::unique_ptr<SolidMechanicsModel> model;
 };
 
 template <typename T> constexpr ElementType TestLumpedMassesFixture<T>::type;
 template <typename T> constexpr size_t TestLumpedMassesFixture<T>::dim;
 
 using mass_types = gtest_list_t<TestElementTypes>;
 
-TYPED_TEST_SUITE(TestLumpedMassesFixture, mass_types);
+TYPED_TEST_SUITE(TestLumpedMassesFixture, mass_types, );
 
 TYPED_TEST(TestLumpedMassesFixture, TestLumpedMass) {
   this->model->assembleMassLumped();
 
   auto rho = this->model->getMaterial(0).getRho();
   auto & fem = this->model->getFEEngine();
   auto nb_element = this->mesh->getNbElement(this->type);
   auto nb_quadrature_points =
       fem.getNbIntegrationPoints(this->type) * nb_element;
 
   Array<Real> rho_on_quad(nb_quadrature_points, 1, rho, "rho_on_quad");
   auto mass = fem.integrate(rho_on_quad, this->type);
   const auto & masses = this->model->getMass();
 
   Vector<Real> sum(this->dim, 0.);
   for (auto & mass : make_view(masses, this->dim)) {
     sum += mass;
   }
 
   for (UInt s = 0; s < sum.size(); ++s)
     EXPECT_NEAR(0., (mass - sum[s]) / mass, 2e-15);
 }
diff --git a/test/test_model/test_common/test_dof_manager.cc b/test/test_model/test_common/test_dof_manager.cc
index a17f2c96a..41cc60eba 100644
--- a/test/test_model/test_common/test_dof_manager.cc
+++ b/test/test_model/test_common/test_dof_manager.cc
@@ -1,298 +1,298 @@
 /**
  * @file   test_dof_manager.cc
  *
  * @author Nicolas Richart
  *
  * @date creation  Wed Jan 30 2019
  *
  * @brief test the dof managers
  *
  *
  * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "test_gtest_utils.hh"
 /* -------------------------------------------------------------------------- */
 #include <dof_manager.hh>
 #include <mesh_partition_scotch.hh>
 #include <mesh_utils.hh>
 /* -------------------------------------------------------------------------- */
 
 #include <gtest/gtest.h>
 #include <numeric>
 #include <string>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 enum DOFManagerType { _dmt_default, _dmt_petsc };
 }
 AKANTU_ENUM_HASH(DOFManagerType)
 
 using namespace akantu;
 
 // defined as struct to get there names in gtest outputs
 struct _dof_manager_default
     : public std::integral_constant<DOFManagerType, _dmt_default> {};
 struct _dof_manager_petsc
     : public std::integral_constant<DOFManagerType, _dmt_petsc> {};
 
 using dof_manager_types = ::testing::Types<
 #ifdef AKANTU_USE_PETSC
     _dof_manager_petsc,
 #endif
     _dof_manager_default>;
 
 namespace std {
 
 std::string to_string(const DOFManagerType & type) {
   std::unordered_map<DOFManagerType, std::string> map{
 #ifdef AKANTU_USE_PETSC
       {_dmt_petsc, "petsc"},
 #endif
       {_dmt_default, "default"},
   };
   return map.at(type);
 }
 
 } // namespace std
 
 /* -------------------------------------------------------------------------- */
 using namespace akantu;
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 class DOFManagerTester {
 public:
   DOFManagerTester(std::unique_ptr<DOFManager> dof_manager)
       : dof_manager(std::move(dof_manager)) {}
 
   DOFManager & operator*() { return *dof_manager; }
   DOFManager * operator->() { return dof_manager.get(); }
   void getArrayPerDOFs(const ID & id, SolverVector & vector,
                        Array<Real> & array) {
     dof_manager->getArrayPerDOFs(id, vector, array);
   }
 
   SolverVector & residual() { return *dof_manager->residual; }
 
 private:
   std::unique_ptr<DOFManager> dof_manager;
 };
 } // namespace akantu
 
 template <class T> class DOFManagerFixture : public ::testing::Test {
 public:
   constexpr static DOFManagerType type = T::value;
   constexpr static UInt dim = 3;
   void SetUp() override {
     mesh = std::make_unique<Mesh>(this->dim);
 
     auto & communicator = Communicator::getStaticCommunicator();
     if (communicator.whoAmI() == 0) {
       mesh->read("mesh.msh");
     }
     mesh->distribute();
 
     nb_nodes = this->mesh->getNbNodes();
     nb_total_nodes = this->mesh->getNbGlobalNodes();
 
     auto && range_nodes = arange(nb_nodes);
     nb_pure_local =
         std::accumulate(range_nodes.begin(), range_nodes.end(), 0,
                         [&](auto && init, auto && val) {
                           return init + mesh->isLocalOrMasterNode(val);
                         });
   }
   void TearDown() override {
     mesh.reset();
     dof1.reset();
     dof2.reset();
   }
 
   decltype(auto) alloc() {
     std::unordered_map<DOFManagerType, std::string> types{
       {_dmt_default, "default"}, {_dmt_petsc, "petsc"}};
 
     return DOFManagerTester(DOFManagerFactory::getInstance().allocate(
         types[T::value], *mesh, "dof_manager", 0));
   }
 
   decltype(auto) registerDOFs(DOFSupportType dst1, DOFSupportType dst2) {
     auto dof_manager = DOFManagerTester(this->alloc());
 
     auto n1 = dst1 == _dst_nodal ? nb_nodes : nb_pure_local;
     this->dof1 = std::make_unique<Array<Real>>(n1, 3);
 
     dof_manager->registerDOFs("dofs1", *this->dof1, dst1);
 
     EXPECT_EQ(dof_manager.residual().size(), nb_total_nodes * 3);
 
     auto n2 = dst2 == _dst_nodal ? nb_nodes : nb_pure_local;
     this->dof2 = std::make_unique<Array<Real>>(n2, 5);
 
     dof_manager->registerDOFs("dofs2", *this->dof2, dst2);
 
     EXPECT_EQ(dof_manager.residual().size(), nb_total_nodes * 8);
     return dof_manager;
   }
 
 protected:
   Int nb_nodes{0}, nb_total_nodes{0}, nb_pure_local{0};
   std::unique_ptr<Mesh> mesh;
   std::unique_ptr<Array<Real>> dof1;
   std::unique_ptr<Array<Real>> dof2;
 };
 
 template <class T> constexpr DOFManagerType DOFManagerFixture<T>::type;
 template <class T> constexpr UInt DOFManagerFixture<T>::dim;
 
-TYPED_TEST_SUITE(DOFManagerFixture, dof_manager_types);
+TYPED_TEST_SUITE(DOFManagerFixture, dof_manager_types, );
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(DOFManagerFixture, Construction) {
   auto dof_manager = this->alloc();
 }
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(DOFManagerFixture, DoubleConstruction) {
   auto dof_manager = this->alloc();
   dof_manager = this->alloc();
 }
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(DOFManagerFixture, RegisterGenericDOF1) {
   auto dof_manager = this->alloc();
 
   Array<Real> dofs(this->nb_pure_local, 3);
 
   dof_manager->registerDOFs("dofs1", dofs, _dst_generic);
   EXPECT_GE(dof_manager.residual().size(), this->nb_total_nodes * 3);
 }
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(DOFManagerFixture, RegisterNodalDOF1) {
   auto dof_manager = this->alloc();
 
   Array<Real> dofs(this->nb_nodes, 3);
   dof_manager->registerDOFs("dofs1", dofs, _dst_nodal);
   EXPECT_GE(dof_manager.residual().size(), this->nb_total_nodes * 3);
 }
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(DOFManagerFixture, RegisterGenericDOF2) {
   this->registerDOFs(_dst_generic, _dst_generic);
 }
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(DOFManagerFixture, RegisterNodalDOF2) {
   this->registerDOFs(_dst_nodal, _dst_nodal);
 }
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(DOFManagerFixture, RegisterMixedDOF) {
   auto dof_manager = this->registerDOFs(_dst_nodal, _dst_generic);
 }
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(DOFManagerFixture, AssembleVector) {
   auto dof_manager = this->registerDOFs(_dst_nodal, _dst_generic);
 
   dof_manager.residual().zero();
 
   for (auto && data :
        enumerate(make_view(*this->dof1, this->dof1->getNbComponent()))) {
     auto n = std::get<0>(data);
     auto & l = std::get<1>(data);
     l.set(1. * this->mesh->isLocalOrMasterNode(n));
   }
 
   this->dof2->set(2.);
 
   dof_manager->assembleToResidual("dofs1", *this->dof1);
   dof_manager->assembleToResidual("dofs2", *this->dof2);
 
   this->dof1->set(0.);
   this->dof2->set(0.);
 
   dof_manager.getArrayPerDOFs("dofs1", dof_manager.residual(), *this->dof1);
   for (auto && data :
        enumerate(make_view(*this->dof1, this->dof1->getNbComponent()))) {
     if (this->mesh->isLocalOrMasterNode(std::get<0>(data))) {
       const auto & l = std::get<1>(data);
       auto e = (l - Vector<Real>{1., 1., 1.}).norm();
       ASSERT_EQ(e, 0.);
     }
   }
 
   dof_manager.getArrayPerDOFs("dofs2", dof_manager.residual(), *this->dof2);
   for (auto && l : make_view(*this->dof2, this->dof2->getNbComponent())) {
     auto e = (l - Vector<Real>{2., 2., 2., 2., 2.}).norm();
     ASSERT_EQ(e, 0.);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(DOFManagerFixture, AssembleMatrixNodal) {
   auto dof_manager = this->registerDOFs(_dst_nodal, _dst_nodal);
 
   auto && K = dof_manager->getNewMatrix("K", _symmetric);
   K.zero();
 
   auto && elemental_matrix = std::make_unique<Array<Real>>(
       this->mesh->getNbElement(this->dim), 8 * 3 * 8 * 3);
 
   for (auto && m : make_view(*elemental_matrix, 8 * 3, 8 * 3)) {
     m.set(1.);
   }
 
   dof_manager->assembleElementalMatricesToMatrix(
       "K", "dofs1", *elemental_matrix, _hexahedron_8);
 
   elemental_matrix = std::make_unique<Array<Real>>(
       this->mesh->getNbElement(this->dim), 8 * 5 * 8 * 5);
 
   for (auto && m : make_view(*elemental_matrix, 8 * 5, 8 * 5)) {
     m.set(1.);
   }
 
   dof_manager->assembleElementalMatricesToMatrix(
       "K", "dofs2", *elemental_matrix, _hexahedron_8);
 
   CSR<Element> node_to_elem;
   MeshUtils::buildNode2Elements(*this->mesh, node_to_elem, this->dim);
 
   dof_manager.residual().zero();
 
   for (auto && data :
        enumerate(zip(make_view(*this->dof1, this->dof1->getNbComponent()),
                      make_view(*this->dof2, this->dof2->getNbComponent())))) {
     auto n = std::get<0>(data);
     auto & l1 = std::get<0>(std::get<1>(data));
     auto & l2 = std::get<1>(std::get<1>(data));
     auto v = 1. * this->mesh->isLocalOrMasterNode(n);
     l1.set(v);
     l2.set(v);
   }
 
   dof_manager->assembleToResidual("dofs1", *this->dof1);
   dof_manager->assembleToResidual("dofs2", *this->dof2);
 
   for (auto && n : arange(this->nb_nodes)) {
     if (not this->mesh->isLocalOrMasterNode(n)) {
     }
   }
 }
diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh
index 596b1655e..d91ddf669 100644
--- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh
+++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh
@@ -1,347 +1,347 @@
 /**
  * @file   test_cohesive_fixture.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Jan 10 2018
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Coehsive element test fixture
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "communicator.hh"
 #include "solid_mechanics_model_cohesive.hh"
 #include "test_gtest_utils.hh"
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_TEST_COHESIVE_FIXTURE_HH_
 #define AKANTU_TEST_COHESIVE_FIXTURE_HH_
 
 using namespace akantu;
 
 template <::akantu::AnalysisMethod t>
 using analysis_method_t = std::integral_constant<::akantu::AnalysisMethod, t>;
 
 class StrainIncrement : public BC::Functor {
 public:
   StrainIncrement(const Matrix<Real> & strain, BC::Axis dir)
       : strain_inc(strain), dir(dir) {}
 
   void operator()(UInt /*node*/, Vector<bool> & flags, Vector<Real> & primal,
                   const Vector<Real> & coord) const {
     if (std::abs(coord(dir)) < 1e-8) {
       return;
     }
 
     flags.set(true);
     primal += strain_inc * coord;
   }
 
   static const BC::Functor::Type type = BC::Functor::_dirichlet;
 
 private:
   Matrix<Real> strain_inc;
   BC::Axis dir;
 };
 
 template <typename param_> class TestSMMCFixture : public ::testing::Test {
 public:
   static constexpr ElementType cohesive_type =
       std::tuple_element_t<0, param_>::value;
   static constexpr ElementType type_1 = std::tuple_element_t<1, param_>::value;
   static constexpr ElementType type_2 = std::tuple_element_t<2, param_>::value;
 
   static constexpr size_t dim =
       ElementClass<cohesive_type>::getSpatialDimension();
 
   void SetUp()  {
     mesh = std::make_unique<Mesh>(this->dim);
     if (Communicator::getStaticCommunicator().whoAmI() == 0) {
        mesh->read(this->mesh_name);
     }
     mesh->distribute();
   }
 
   void TearDown() {
     model.reset(nullptr);
     mesh.reset(nullptr);
   }
 
   void createModel() {
     model = std::make_unique<SolidMechanicsModelCohesive>(*mesh);
     model->initFull(_analysis_method = this->analysis_method,
                     _is_extrinsic = this->is_extrinsic);
 
     auto time_step = this->model->getStableTimeStep() * 0.01;
     this->model->setTimeStep(time_step);
 
     if (dim == 1) {
       surface = 1;
       group_size = 1;
       return;
     }
 
     auto facet_type = mesh->getFacetType(this->cohesive_type);
 
     auto & fe_engine = model->getFEEngineBoundary();
     const auto & group = mesh->getElementGroup("insertion");
     group_size = group.size(_ghost_type = _not_ghost);
     const auto & elements = group.getElements(facet_type);
     Array<Real> ones(fe_engine.getNbIntegrationPoints(facet_type) *
                      group_size);
     ones.set(1.);
 
     surface = fe_engine.integrate(ones, facet_type, _not_ghost, elements);
     mesh->getCommunicator().allReduce(surface, SynchronizerOperation::_sum);
 
 
 
     mesh->getCommunicator().allReduce(group_size, SynchronizerOperation::_sum);
 
 #define debug_ 0
 
 #if debug_
     this->model->addDumpFieldVector("displacement");
     this->model->addDumpFieldVector("velocity");
     this->model->addDumpFieldVector("internal_force");
     this->model->addDumpFieldVector("external_force");
     this->model->addDumpField("blocked_dofs");
     this->model->addDumpField("stress");
     this->model->addDumpField("strain");
     this->model->assembleInternalForces();
     this->model->setBaseNameToDumper("cohesive elements", "cohesive_elements");
     this->model->addDumpFieldVectorToDumper("cohesive elements",
                                             "displacement");
     this->model->addDumpFieldToDumper("cohesive elements", "damage");
     this->model->addDumpFieldToDumper("cohesive elements", "tractions");
     this->model->addDumpFieldToDumper("cohesive elements", "opening");
     this->model->dump();
     this->model->dump("cohesive elements");
 #endif
   }
 
   void setInitialCondition(const Matrix<Real> & strain) {
     for (auto && data :
          zip(make_view(this->mesh->getNodes(), this->dim),
              make_view(this->model->getDisplacement(), this->dim))) {
       const auto & pos = std::get<0>(data);
       auto & disp = std::get<1>(data);
       disp = strain * pos;
     }
   }
 
   bool checkDamaged() {
     UInt nb_damaged = 0;
     auto & damage =
         model->getMaterial("insertion").getArray<Real>("damage", cohesive_type);
     for (auto d : damage) {
       if (d >= .99) {
         ++nb_damaged;
       }
     }
 
     return (nb_damaged == group_size);
   }
 
   void steps(const Matrix<Real> & strain) {
     StrainIncrement functor((1. / 300) * strain, this->dim == 1 ? _x : _y);
 
     for (auto _ [[gnu::unused]] : arange(nb_steps)) {
       this->model->applyBC(functor, "loading");
       this->model->applyBC(functor, "fixed");
       if (this->is_extrinsic) {
         this->model->checkCohesiveStress();
       }
 
       this->model->solveStep();
 #if debug_
       this->model->dump();
       this->model->dump("cohesive elements");
 #endif
     }
   }
 
   void checkInsertion() {
     auto nb_cohesive_element = this->mesh->getNbElement(cohesive_type);
     mesh->getCommunicator().allReduce(nb_cohesive_element,
                                       SynchronizerOperation::_sum);
 
     EXPECT_EQ(nb_cohesive_element, group_size);
   }
 
   void checkDissipated(Real expected_density) {
     Real edis = this->model->getEnergy("dissipated");
 
     EXPECT_NEAR(this->surface * expected_density, edis, 5e-1);
   }
 
   void testModeI() {
     this->createModel();
 
     auto & mat_el = this->model->getMaterial("body");
 
     auto speed = mat_el.getPushWaveSpeed(Element());
     auto direction = _y;
     if (dim == 1) {
       direction = _x;
     }
     auto length =
         mesh->getUpperBounds()(direction) - mesh->getLowerBounds()(direction);
     nb_steps = length / speed / model->getTimeStep();
 
     SCOPED_TRACE(std::to_string(this->dim) + "D - " + std::to_string(type_1) +
                  ":" + std::to_string(type_2));
 
     auto & mat_co = this->model->getMaterial("insertion");
     Real sigma_c = mat_co.get("sigma_c");
 
     Real E = mat_el.get("E");
     Real nu = mat_el.get("nu");
 
     Matrix<Real> strain;
     if (dim == 1) {
       strain = {{1.}};
     } else if (dim == 2) {
       strain = {{-nu, 0.}, {0., 1. - nu}};
       strain *= (1. + nu);
     } else if (dim == 3) {
       strain = {{-nu, 0., 0.}, {0., 1., 0.}, {0., 0., -nu}};
     }
 
     strain *= sigma_c / E;
 
     this->setInitialCondition((1 - 1e-5) * strain);
     this->steps(2e-2 * strain);
   }
 
   void testModeII() {
     this->createModel();
     auto & mat_el = this->model->getMaterial("body");
     Real speed;
     try {
       speed =
           mat_el.getShearWaveSpeed(Element()); // the slowest speed if exists
     } catch (...) {
       speed = mat_el.getPushWaveSpeed(Element());
     }
 
     auto direction = _y;
     if (dim == 1)
       direction = _x;
     auto length =
         mesh->getUpperBounds()(direction) - mesh->getLowerBounds()(direction);
     nb_steps = 2 * length / 2. / speed / model->getTimeStep();
 
     SCOPED_TRACE(std::to_string(this->dim) + "D - " + std::to_string(type_1) +
                  ":" + std::to_string(type_2));
 
     if (this->dim > 1)
       this->model->applyBC(BC::Dirichlet::FlagOnly(_y), "sides");
     if (this->dim > 2)
       this->model->applyBC(BC::Dirichlet::FlagOnly(_z), "sides");
 
     auto & mat_co = this->model->getMaterial("insertion");
     Real sigma_c = mat_co.get("sigma_c");
     Real beta = mat_co.get("beta");
     // Real G_c = mat_co.get("G_c");
 
     Real E = mat_el.get("E");
     Real nu = mat_el.get("nu");
 
     Matrix<Real> strain;
     if (dim == 1) {
       strain = {{1.}};
     } else if (dim == 2) {
       strain = {{0., 1.}, {0., 0.}};
       strain *= (1. + nu);
     } else if (dim == 3) {
       strain = {{0., 1., 0.}, {0., 0., 0.}, {0., 0., 0.}};
       strain *= (1. + nu);
     }
     strain *= 2 * beta * beta * sigma_c / E;
 
     // nb_steps *= 5;
 
     this->setInitialCondition((1. - 1e-5) * strain);
     this->steps(0.005 * strain);
   }
 
 protected:
   std::unique_ptr<Mesh> mesh;
   std::unique_ptr<SolidMechanicsModelCohesive> model;
 
   std::string mesh_name{std::to_string(cohesive_type) + std::to_string(type_1) +
                         (type_1 == type_2 ? "" : std::to_string(type_2)) +
                         ".msh"};
 
   bool is_extrinsic;
   AnalysisMethod analysis_method;
 
   Real surface{0};
   UInt nb_steps{1000};
   UInt group_size{10000};
 };
 
 /* -------------------------------------------------------------------------- */
 template <typename param_>
 constexpr ElementType TestSMMCFixture<param_>::cohesive_type;
 template <typename param_>
 constexpr ElementType TestSMMCFixture<param_>::type_1;
 template <typename param_>
 constexpr ElementType TestSMMCFixture<param_>::type_2;
 
 template <typename param_> constexpr size_t TestSMMCFixture<param_>::dim;
 /* -------------------------------------------------------------------------- */
 
 using IsExtrinsicTypes = std::tuple<std::true_type, std::false_type>;
 using AnalysisMethodTypes =
     std::tuple<analysis_method_t<_explicit_lumped_mass>>;
 
 using coh_types = gtest_list_t<std::tuple<
     std::tuple<_element_type_cohesive_1d_2, _element_type_segment_2,
                _element_type_segment_2>,
     std::tuple<_element_type_cohesive_2d_4, _element_type_triangle_3,
                _element_type_triangle_3>,
     std::tuple<_element_type_cohesive_2d_4, _element_type_quadrangle_4,
                _element_type_quadrangle_4>,
     std::tuple<_element_type_cohesive_2d_4, _element_type_triangle_3,
                _element_type_quadrangle_4>,
     std::tuple<_element_type_cohesive_2d_6, _element_type_triangle_6,
                _element_type_triangle_6>,
     std::tuple<_element_type_cohesive_2d_6, _element_type_quadrangle_8,
                _element_type_quadrangle_8>,
     std::tuple<_element_type_cohesive_2d_6, _element_type_triangle_6,
                _element_type_quadrangle_8>,
     std::tuple<_element_type_cohesive_3d_6, _element_type_tetrahedron_4,
                _element_type_tetrahedron_4>,
     std::tuple<_element_type_cohesive_3d_12, _element_type_tetrahedron_10,
                _element_type_tetrahedron_10> /*,
     std::tuple<_element_type_cohesive_3d_8, _element_type_hexahedron_8,
                _element_type_hexahedron_8>,
     std::tuple<_element_type_cohesive_3d_16, _element_type_hexahedron_20,
                _element_type_hexahedron_20>*/>>;
 
-TYPED_TEST_SUITE(TestSMMCFixture, coh_types);
+TYPED_TEST_SUITE(TestSMMCFixture, coh_types, );
 
 #endif /* AKANTU_TEST_COHESIVE_FIXTURE_HH_ */
diff --git a/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc b/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc
index 102020166..95bb94ec0 100644
--- a/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc
+++ b/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc
@@ -1,157 +1,157 @@
 /**
  * @file   test_solid_mechanics_model_work_dynamics.cc
  *
  * @author Tobias Brink <tobias.brink@epfl.ch>
  *
  * @date creation: Fri Dec 15 2017
  * @date last modification: Fri Jan 26 2018
  *
  * @brief  test work in dynamic simulations
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  * @section description
  *
  * Assuming that the kinetic energy and the potential energy of the
  * linear elastic material are bug free, the work in a dynamic
  * simulation must equal the change in internal energy (first law of
  * thermodynamics). Work in dynamics is an infinitesimal work Fds,
  * thus we need to integrate it and compare at the end. In this test,
  * we use one Dirichlet boundary condition (with u = 0.0, 0.01, and
  * -0.01) and one Neumann boundary condition for F on the opposite
  * side. Then we do a few steps to get reference energies for work and
  * internal energy. After more steps, the change in both work and
  * internal energy must be equal.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "../test_solid_mechanics_model_fixture.hh"
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 namespace {
 
 template <typename type_>
 class TestSMMFixtureWorkDynamic : public TestSMMFixture<type_> {
 public:
   void SetUp() override {
     this->mesh_file =
         "../../patch_tests/data/bar" + std::to_string(this->type) + ".msh";
     TestSMMFixture<type_>::SetUp();
 
     getStaticParser().parse("test_solid_mechanics_model_"
                             "work_material.dat");
 
     /// model initialization
     this->model->initFull();
 
     /// Create a node group for Neumann BCs.
     auto & apply_force_grp = this->mesh->createNodeGroup("apply_force");
     auto & fixed_grp = this->mesh->createNodeGroup("fixed");
 
     const auto & pos = this->mesh->getNodes();
     const auto & lower = this->mesh->getLowerBounds();
     const auto & upper = this->mesh->getUpperBounds();
 
     UInt i = 0;
     for (auto && posv : make_view(pos, this->spatial_dimension)) {
       if (posv(_x) > upper(_x) - 1e-6) {
         apply_force_grp.add(i);
       } else if (posv(_x) < lower(_x) + 1e-6) {
         fixed_grp.add(i);
       }
       ++i;
     }
 
     this->mesh->createElementGroupFromNodeGroup("el_apply_force", "apply_force",
                                                 this->spatial_dimension - 1);
     this->mesh->createElementGroupFromNodeGroup("el_fixed", "fixed",
                                                 this->spatial_dimension - 1);
 
     Vector<Real> surface_traction(this->spatial_dimension);
     surface_traction(_x) = 0.5;
     if (this->spatial_dimension == 1) {
       // TODO: this is a hack to work
       //      around non-implemented
       //      BC::Neumann::FromTraction for 1D
       auto & force = this->model->getExternalForce();
 
       for (auto && pair : zip(make_view(pos, this->spatial_dimension),
                               make_view(force, this->spatial_dimension))) {
         auto & posv = std::get<0>(pair);
         auto & forcev = std::get<1>(pair);
         if (posv(_x) > upper(_x) - 1e-6) {
           forcev(_x) = surface_traction(_x);
         }
       }
     } else {
       this->model->applyBC(BC::Neumann::FromTraction(surface_traction),
                            "el_apply_force");
     }
 
     /// set up timestep
     auto time_step = this->model->getStableTimeStep() * 0.1;
     this->model->setTimeStep(time_step);
   }
 };
 
-TYPED_TEST_SUITE(TestSMMFixtureWorkDynamic, gtest_element_types);
+TYPED_TEST_SUITE(TestSMMFixtureWorkDynamic, gtest_element_types, );
 
 /* TODO: this is currently disabled for terrible results and performance
 TYPED_TEST(TestSMMFixtureBar, WorkImplicit) {
   test_body(*(this->model), *(this->mesh), _implicit_dynamic, 500);
 }
 */
 // model.assembleMassLumped();
 TYPED_TEST(TestSMMFixtureWorkDynamic, WorkExplicit) {
   /// Do the sim
   std::vector<Real> displacements{0.00, 0.01, -0.01};
 
   for (auto && u : displacements) {
     this->model->applyBC(BC::Dirichlet::FixedValue(u, _x), "el_fixed");
 
     // First, "equilibrate" a bit to get a reference state of total
     // energy and work. This is needed when we have a Dirichlet with
     // finite displacement on one side.
     for (UInt i = 0; i < 25; ++i) {
       this->model->solveStep();
     }
     // Again, work reported by Akantu is infinitesimal (dW) and we
     // need to integrate a while to get a decent value.
     double Etot0 =
         this->model->getEnergy("potential") + this->model->getEnergy("kinetic");
     double W = 0.0;
     for (UInt i = 0; i < 200; ++i) {
       /// Solve.
       this->model->solveStep();
 
       const auto dW = this->model->getEnergy("external work");
 
       W += dW;
     }
     // Finally check.
     const auto Epot = this->model->getEnergy("potential");
     const auto Ekin = this->model->getEnergy("kinetic");
 
     EXPECT_NEAR(W, Ekin + Epot - Etot0, 5e-2);
     // Sadly not very exact for such a coarse mesh.
   }
 }
 } // namespace
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_damage_materials.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_damage_materials.cc
index 82b0d4b3e..7e3e9c4ac 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_damage_materials.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_damage_materials.cc
@@ -1,243 +1,243 @@
 /**
  * @file   test_damage_materials.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @date creation: Fri Nov 17 2017
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Tests for damage materials
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "test_material_fixtures.hh"
 /* -------------------------------------------------------------------------- */
 #include <material_marigo.hh>
 #include <material_mazars.hh>
 #include <py_aka_array.hh>
 #include <solid_mechanics_model.hh>
 /* -------------------------------------------------------------------------- */
 #include <fstream>
 #include <gtest/gtest.h>
 #include <pybind11/embed.h>
 #include <pybind11/numpy.h>
 #include <pybind11/stl.h>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 namespace py = pybind11;
 using namespace py::literals;
 
 using mat_types = ::testing::Types<
     // Traits<MaterialMarigo, 1>, Traits<MaterialMarigo, 2>,
     // Traits<MaterialMarigo, 3>,
     Traits<MaterialMazars, 1>, Traits<MaterialMazars, 2>,
     Traits<MaterialMazars, 3>>;
 
 /*****************************************************************/
 
 template <> void FriendMaterial<MaterialMazars<1>>::setParams() {
   K0.setDefaultValue(1e-4);
   At = 1.0;
   Bt = 5e3;
   Ac = 0.8;
   Bc = 1391.3;
   beta = 1.;
   E = 25e9;
   nu = 0.2;
 
   updateInternalParameters();
 }
 
 template <> void FriendMaterial<MaterialMazars<2>>::setParams() {
   K0.setDefaultValue(1e-4);
   At = 1.0;
   Bt = 5e3;
   Ac = 0.8;
   Bc = 1391.3;
   beta = 1.;
   E = 25e9;
   nu = 0.2;
   plane_stress = true;
   updateInternalParameters();
 }
 
 template <> void FriendMaterial<MaterialMazars<3>>::setParams() {
   K0.setDefaultValue(1e-4);
   At = 1.0;
   Bt = 5e3;
   Ac = 0.8;
   Bc = 1391.3;
   beta = 1.;
   E = 25e9;
   nu = 0.2;
 
   updateInternalParameters();
 }
 
 template <> void FriendMaterial<MaterialMazars<1>>::testComputeStress() {
   Array<Real> epsilons(1001, 1);
   Array<Real> sigmas(1001, 1);
   Array<Real> damages(1001, 1);
 
   for (auto && data : enumerate(epsilons)) {
     std::get<1>(data) = 2e-6 * std::get<0>(data);
   }
   Real _K0 = K0;
   py::module py_engine = py::module::import("py_mazars");
 
   auto kwargs_mat_params =
       py::dict("K0"_a = _K0, "At"_a = At, "Bt"_a = Bt, "Ac"_a = Ac, "Bc"_a = Bc,
                "E"_a = E, "nu"_a = nu);
   auto kwargs = py::dict("epsilons"_a = epsilons, "sigmas"_a = sigmas,
                          "damages"_a = damages);
 
   auto py_mazars = py_engine.attr("Mazars")(**kwargs_mat_params);
   // auto Gf_py = py_mazars.attr("compute")(**kwargs);
 
   Real dam = 0.;
   Real dam_ref = 0.;
   Real ehat = 0.;
 
   for (auto && epsilon : epsilons) {
     Matrix<Real> strain(this->spatial_dimension, this->spatial_dimension, 0.);
     Matrix<Real> sigma(this->spatial_dimension, this->spatial_dimension, 0.);
     strain(0, 0) = epsilon;
 
     computeStressOnQuad(strain, sigma, dam, ehat);
 
     Real sigma_ref;
     auto py_data =
         py_mazars.attr("compute_step")(epsilon, sigma_ref, dam_ref, false);
     std::tie(sigma_ref, dam_ref) = py::cast<std::pair<double, double>>(py_data);
 
     EXPECT_NEAR(sigma(0, 0), sigma_ref, 1e-5);
     EXPECT_NEAR(dam, dam_ref, 1e-10);
   }
 }
 
 template <> void FriendMaterial<MaterialMazars<2>>::testComputeStress() {
   Array<Real> epsilons(1001, 1);
   Array<Real> sigmas(1001, 1);
   Array<Real> damages(1001, 1);
 
   for (auto && data : enumerate(epsilons)) {
     std::get<1>(data) = 2e-6 * std::get<0>(data);
   }
   Real _K0 = K0;
   py::module py_engine = py::module::import("py_mazars");
 
   auto kwargs_mat_params =
       py::dict("K0"_a = _K0, "At"_a = At, "Bt"_a = Bt, "Ac"_a = Ac, "Bc"_a = Bc,
                "E"_a = E, "nu"_a = nu);
   auto kwargs = py::dict("epsilons"_a = epsilons, "sigmas"_a = sigmas,
                          "damages"_a = damages);
 
   auto py_mazars = py_engine.attr("Mazars")(**kwargs_mat_params);
   // auto Gf_py = py_mazars.attr("compute")(**kwargs);
 
   Real dam = 0.;
   Real dam_ref = 0.;
   Real ehat = 0.;
 
   for (auto && epsilon : epsilons) {
     Matrix<Real> strain(this->spatial_dimension, this->spatial_dimension, 0.);
     Matrix<Real> sigma(this->spatial_dimension, this->spatial_dimension, 0.);
     strain(0, 0) = epsilon;
     strain(1, 1) = -this->nu * epsilon;
 
     computeStressOnQuad(strain, sigma, dam, ehat);
 
     Real sigma_ref;
     auto py_data =
         py_mazars.attr("compute_step")(epsilon, sigma_ref, dam_ref, false);
     std::tie(sigma_ref, dam_ref) = py::cast<std::pair<double, double>>(py_data);
 
     EXPECT_NEAR(sigma(0, 0), sigma_ref, 1e-5);
     EXPECT_NEAR(dam, dam_ref, 1e-10);
   }
 }
 
 template <> void FriendMaterial<MaterialMazars<3>>::testComputeStress() {
   Array<Real> epsilons(1001, 1);
   Array<Real> sigmas(1001, 1);
   Array<Real> damages(1001, 1);
 
   for (auto && data : enumerate(epsilons)) {
     std::get<1>(data) = 2e-6 * std::get<0>(data);
   }
   Real _K0 = K0;
   py::module py_engine = py::module::import("py_mazars");
 
   auto kwargs_mat_params =
       py::dict("K0"_a = _K0, "At"_a = At, "Bt"_a = Bt, "Ac"_a = Ac, "Bc"_a = Bc,
                "E"_a = E, "nu"_a = nu);
   auto kwargs = py::dict("epsilons"_a = epsilons, "sigmas"_a = sigmas,
                          "damages"_a = damages);
 
   auto py_mazars = py_engine.attr("Mazars")(**kwargs_mat_params);
   // auto Gf_py = py_mazars.attr("compute")(**kwargs);
 
   Real dam = 0.;
   Real dam_ref = 0.;
   Real ehat = 0.;
 
   for (auto && epsilon : epsilons) {
     Matrix<Real> strain(this->spatial_dimension, this->spatial_dimension, 0.);
     Matrix<Real> sigma(this->spatial_dimension, this->spatial_dimension, 0.);
     strain(0, 0) = epsilon;
     strain(1, 1) = strain(2, 2) = -this->nu * epsilon;
 
     computeStressOnQuad(strain, sigma, dam, ehat);
 
     Real sigma_ref;
     auto py_data =
         py_mazars.attr("compute_step")(epsilon, sigma_ref, dam_ref, false);
     std::tie(sigma_ref, dam_ref) = py::cast<std::pair<double, double>>(py_data);
 
     EXPECT_NEAR(sigma(0, 0), sigma_ref, 1e-5);
     EXPECT_NEAR(dam, dam_ref, 1e-10);
   }
 }
 
 namespace {
 
 template <typename T>
 class TestDamageMaterialFixture : public ::TestMaterialFixture<T> {};
 
-TYPED_TEST_SUITE(TestDamageMaterialFixture, mat_types);
+TYPED_TEST_SUITE(TestDamageMaterialFixture, mat_types, );
 
 TYPED_TEST(TestDamageMaterialFixture, ComputeStress) {
   this->material->testComputeStress();
 }
 TYPED_TEST(TestDamageMaterialFixture, DISABLED_EnergyDensity) {
   this->material->testEnergyDensity();
 }
 TYPED_TEST(TestDamageMaterialFixture, DISABLED_ComputeTangentModuli) {
   this->material->testComputeTangentModuli();
 }
 TYPED_TEST(TestDamageMaterialFixture, DISABLED_ComputeCelerity) {
   this->material->testCelerity();
 }
 } // namespace
 /*****************************************************************/
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc
index ff5f308a5..9e4d05a61 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc
@@ -1,884 +1,884 @@
 /**
  * @file   test_elastic_materials.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  * @author Enrico Milanese <enrico.milanese@epfl.ch>
  *
  * @date creation: Fri Nov 17 2017
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Tests the Elastic materials
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "test_gtest_utils.hh"
 #include "test_material_fixtures.hh"
 /* -------------------------------------------------------------------------- */
 #include <material_elastic.hh>
 #include <material_elastic_orthotropic.hh>
 #include <solid_mechanics_model.hh>
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 using mat_types =
     ::testing::Types<Traits<MaterialElastic, 1>, Traits<MaterialElastic, 2>,
                      Traits<MaterialElastic, 3>,
 
                      Traits<MaterialElasticOrthotropic, 2>,
                      Traits<MaterialElasticOrthotropic, 3>,
 
                      Traits<MaterialElasticLinearAnisotropic, 2>,
                      Traits<MaterialElasticLinearAnisotropic, 3>>;
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<1>>::setParams() {
   Real E = 3.;
   Real rho = 2;
   setParam("E", E);
   setParam("rho", rho);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<1>>::testComputeStress() {
   Matrix<Real> eps = {{2}};
   Matrix<Real> sigma(1, 1);
   Real sigma_th = 2;
   this->computeStressOnQuad(eps, sigma, sigma_th);
 
   auto solution = E * eps(0, 0) + sigma_th;
   EXPECT_NEAR(sigma(0, 0), solution, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<1>>::testEnergyDensity() {
   Real eps = 2, sigma = 2;
   Real epot = 0;
   this->computePotentialEnergyOnQuad({{eps}}, {{sigma}}, epot);
   Real solution = 2;
   EXPECT_NEAR(epot, solution, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElastic<1>>::testComputeTangentModuli() {
   Matrix<Real> tangent(1, 1);
   this->computeTangentModuliOnQuad(tangent);
   EXPECT_NEAR(tangent(0, 0), E, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<1>>::testCelerity() {
   auto wave_speed = this->getCelerity(Element());
   auto solution = std::sqrt(E / rho);
   EXPECT_NEAR(wave_speed, solution, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<2>>::setParams() {
   Real E = 1.;
   Real nu = .3;
   Real rho = 2;
   setParam("E", E);
   setParam("nu", nu);
   setParam("rho", rho);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<2>>::testComputeStress() {
   Real bulk_modulus_K = E / (3 * (1 - 2 * nu));
   Real shear_modulus_mu = E / (2 * (1 + nu));
 
   auto rotation_matrix = getRandomRotation();
 
   auto grad_u = this->getComposedStrain(1.).block(0, 0, 2, 2);
 
   auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix);
 
   Matrix<Real> sigma_rot(2, 2);
   this->computeStressOnQuad(grad_u_rot, sigma_rot, sigma_th);
 
   auto sigma = this->reverseRotation(sigma_rot, rotation_matrix);
 
   auto identity = Matrix<Real>::eye(2, 1.);
 
   auto strain = 0.5 * (grad_u + grad_u.transpose());
   auto deviatoric_strain = strain - 1. / 3. * strain.trace() * identity;
 
   auto sigma_expected = 2 * shear_modulus_mu * deviatoric_strain +
                         (sigma_th + 2. * bulk_modulus_K) * identity;
 
   auto diff = sigma - sigma_expected;
   Real stress_error = diff.norm<L_inf>() / sigma_expected.norm<L_inf>();
   EXPECT_NEAR(stress_error, 0., 1e-13);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<2>>::testEnergyDensity() {
   Matrix<Real> sigma = {{1, 2}, {2, 4}};
   Matrix<Real> eps = {{1, 0}, {0, 1}};
   Real epot = 0;
   Real solution = 2.5;
   this->computePotentialEnergyOnQuad(eps, sigma, epot);
   EXPECT_NEAR(epot, solution, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElastic<2>>::testComputeTangentModuli() {
   Matrix<Real> tangent(3, 3);
 
   /* Plane Strain */
   // clang-format off
   Matrix<Real> solution = {
     {1 - nu, nu, 0},
     {nu, 1 - nu, 0},
     {0, 0, (1 - 2 * nu) / 2},
   };
   // clang-format on
   solution *= E / ((1 + nu) * (1 - 2 * nu));
 
   this->computeTangentModuliOnQuad(tangent);
   Real tangent_error = (tangent - solution).norm<L_2>();
   EXPECT_NEAR(tangent_error, 0, 1e-14);
 
   /* Plane Stress */
   this->plane_stress = true;
   this->updateInternalParameters();
   // clang-format off
   solution = {
     {1, nu, 0},
     {nu, 1, 0},
     {0, 0, (1 - nu) / 2},
   };
   // clang-format on
   solution *= E / (1 - nu * nu);
 
   this->computeTangentModuliOnQuad(tangent);
   tangent_error = (tangent - solution).norm<L_2>();
   EXPECT_NEAR(tangent_error, 0, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<2>>::testCelerity() {
   auto push_wave_speed = this->getPushWaveSpeed(Element());
   auto celerity = this->getCelerity(Element());
 
   Real K = E / (3 * (1 - 2 * nu));
   Real mu = E / (2 * (1 + nu));
   Real sol = std::sqrt((K + 4. / 3 * mu) / rho);
 
   EXPECT_NEAR(push_wave_speed, sol, 1e-14);
   EXPECT_NEAR(celerity, sol, 1e-14);
 
   auto shear_wave_speed = this->getShearWaveSpeed(Element());
 
   sol = std::sqrt(mu / rho);
 
   EXPECT_NEAR(shear_wave_speed, sol, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<3>>::setParams() {
   Real E = 1.;
   Real nu = .3;
   Real rho = 2;
   setParam("E", E);
   setParam("nu", nu);
   setParam("rho", rho);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<3>>::testComputeStress() {
   Real bulk_modulus_K = E / 3. / (1 - 2. * nu);
   Real shear_modulus_mu = 0.5 * E / (1 + nu);
 
   Matrix<Real> rotation_matrix = getRandomRotation();
 
   auto grad_u = this->getComposedStrain(1.);
 
   auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix);
 
   Matrix<Real> sigma_rot(3, 3);
   this->computeStressOnQuad(grad_u_rot, sigma_rot, sigma_th);
 
   auto sigma = this->reverseRotation(sigma_rot, rotation_matrix);
 
   Matrix<Real> identity(3, 3);
   identity.eye();
 
   Matrix<Real> strain = 0.5 * (grad_u + grad_u.transpose());
   Matrix<Real> deviatoric_strain = strain - 1. / 3. * strain.trace() * identity;
 
   Matrix<Real> sigma_expected = 2 * shear_modulus_mu * deviatoric_strain +
                                 (sigma_th + 3. * bulk_modulus_K) * identity;
 
   auto diff = sigma - sigma_expected;
   Real stress_error = diff.norm<L_inf>();
   EXPECT_NEAR(stress_error, 0., 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<3>>::testEnergyDensity() {
   Matrix<Real> sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
   Matrix<Real> eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
   Real epot = 0;
   Real solution = 5.5;
   this->computePotentialEnergyOnQuad(eps, sigma, epot);
   EXPECT_NEAR(epot, solution, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElastic<3>>::testComputeTangentModuli() {
   Matrix<Real> tangent(6, 6);
 
   // clang-format off
   Matrix<Real> solution = {
       {1 - nu, nu, nu, 0, 0, 0},
       {nu, 1 - nu, nu, 0, 0, 0},
       {nu, nu, 1 - nu, 0, 0, 0},
       {0, 0, 0, (1 - 2 * nu) / 2, 0, 0},
       {0, 0, 0, 0, (1 - 2 * nu) / 2, 0},
       {0, 0, 0, 0, 0, (1 - 2 * nu) / 2},
   };
   // clang-format on
   solution *= E / ((1 + nu) * (1 - 2 * nu));
 
   this->computeTangentModuliOnQuad(tangent);
   Real tangent_error = (tangent - solution).norm<L_2>();
   EXPECT_NEAR(tangent_error, 0, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<3>>::testCelerity() {
   auto push_wave_speed = this->getPushWaveSpeed(Element());
   auto celerity = this->getCelerity(Element());
 
   Real K = E / (3 * (1 - 2 * nu));
   Real mu = E / (2 * (1 + nu));
   Real sol = std::sqrt((K + 4. / 3 * mu) / rho);
 
   EXPECT_NEAR(push_wave_speed, sol, 1e-14);
   EXPECT_NEAR(celerity, sol, 1e-14);
 
   auto shear_wave_speed = this->getShearWaveSpeed(Element());
 
   sol = std::sqrt(mu / rho);
 
   EXPECT_NEAR(shear_wave_speed, sol, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::setParams() {
   // Note: for this test material and canonical basis coincide
   Vector<Real> n1 = {1, 0};
   Vector<Real> n2 = {0, 1};
   Real E1 = 1.;
   Real E2 = 2.;
   Real nu12 = 0.1;
   Real G12 = 2.;
   Real rho = 2.5;
 
   *this->dir_vecs[0] = n1;
   *this->dir_vecs[1] = n2;
   this->E1 = E1;
   this->E2 = E2;
   this->nu12 = nu12;
   this->G12 = G12;
   this->rho = rho;
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticOrthotropic<2>>::testComputeStress() {
   UInt Dim = 2;
   // material frame of reference is rotate by rotation_matrix starting from
   // canonical basis
   Matrix<Real> rotation_matrix = getRandomRotation();
 
   // canonical basis as expressed in the material frame of reference, as
   // required by MaterialElasticOrthotropic class (it is simply given by the
   // columns of the rotation_matrix; the lines give the material basis expressed
   // in the canonical frame of reference)
   *this->dir_vecs[0] = rotation_matrix(0);
   *this->dir_vecs[1] = rotation_matrix(1);
 
   // set internal Cijkl matrix expressed in the canonical frame of reference
   this->updateInternalParameters();
 
   // gradient in material frame of reference
   auto grad_u = this->getComposedStrain(2.).block(0, 0, 2, 2);
 
   // gradient in canonical basis (we need to rotate *back* to the canonical
   // basis)
   auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);
 
   // stress in the canonical basis
   Matrix<Real> sigma_rot(2, 2);
   this->computeStressOnQuad(grad_u_rot, sigma_rot);
 
   // stress in the material reference (we need to apply the rotation)
   auto sigma = this->applyRotation(sigma_rot, rotation_matrix);
 
   // construction of Cijkl engineering tensor in the *material* frame of
   // reference
   // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
   Real nu21 = nu12 * E2 / E1;
   Real gamma = 1 / (1 - nu12 * nu21);
 
   Matrix<Real> C_expected(2 * Dim, 2 * Dim, 0);
   C_expected(0, 0) = gamma * E1;
   C_expected(1, 1) = gamma * E2;
   C_expected(2, 2) = G12;
 
   C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * nu21;
 
   // epsilon is computed directly in the *material* frame of reference
   Matrix<Real> epsilon = 0.5 * (grad_u + grad_u.transpose());
 
   // sigma_expected is computed directly in the *material* frame of reference
   Matrix<Real> sigma_expected(Dim, Dim);
   for (UInt i = 0; i < Dim; ++i) {
     for (UInt j = 0; j < Dim; ++j) {
       sigma_expected(i, i) += C_expected(i, j) * epsilon(j, j);
     }
   }
 
   sigma_expected(0, 1) = sigma_expected(1, 0) =
       C_expected(2, 2) * 2 * epsilon(0, 1);
 
   // sigmas are checked in the *material* frame of reference
   auto diff = sigma - sigma_expected;
   Real stress_error = diff.norm<L_inf>();
   EXPECT_NEAR(stress_error, 0., 1e-13);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticOrthotropic<2>>::testEnergyDensity() {
   Matrix<Real> sigma = {{1, 2}, {2, 4}};
   Matrix<Real> eps = {{1, 0}, {0, 1}};
   Real epot = 0;
   Real solution = 2.5;
   this->computePotentialEnergyOnQuad(eps, sigma, epot);
   EXPECT_NEAR(epot, solution, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticOrthotropic<2>>::testComputeTangentModuli() {
   // construction of Cijkl engineering tensor in the *material* frame of
   // reference
   // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
   Real nu21 = nu12 * E2 / E1;
   Real gamma = 1 / (1 - nu12 * nu21);
 
   Matrix<Real> C_expected(3, 3);
   C_expected(0, 0) = gamma * E1;
   C_expected(1, 1) = gamma * E2;
   C_expected(2, 2) = G12;
 
   C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * nu21;
 
   Matrix<Real> tangent(3, 3);
   this->computeTangentModuliOnQuad(tangent);
 
   Real tangent_error = (tangent - C_expected).norm<L_2>();
   EXPECT_NEAR(tangent_error, 0, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testCelerity() {
   // construction of Cijkl engineering tensor in the *material* frame of
   // reference
   // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
   Real nu21 = nu12 * E2 / E1;
   Real gamma = 1 / (1 - nu12 * nu21);
 
   Matrix<Real> C_expected(3, 3);
   C_expected(0, 0) = gamma * E1;
   C_expected(1, 1) = gamma * E2;
   C_expected(2, 2) = G12;
 
   C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * nu21;
 
   Vector<Real> eig_expected(3);
   C_expected.eig(eig_expected);
 
   auto celerity_expected = std::sqrt(eig_expected(0) / rho);
 
   auto celerity = this->getCelerity(Element());
 
   EXPECT_NEAR(celerity_expected, celerity, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::setParams() {
   Vector<Real> n1 = {1, 0, 0};
   Vector<Real> n2 = {0, 1, 0};
   Vector<Real> n3 = {0, 0, 1};
   Real E1 = 1.;
   Real E2 = 2.;
   Real E3 = 3.;
   Real nu12 = 0.1;
   Real nu13 = 0.2;
   Real nu23 = 0.3;
   Real G12 = 2.;
   Real G13 = 3.;
   Real G23 = 1.;
   Real rho = 2.3;
 
   *this->dir_vecs[0] = n1;
   *this->dir_vecs[1] = n2;
   *this->dir_vecs[2] = n3;
   this->E1 = E1;
   this->E2 = E2;
   this->E3 = E3;
   this->nu12 = nu12;
   this->nu13 = nu13;
   this->nu23 = nu23;
   this->G12 = G12;
   this->G13 = G13;
   this->G23 = G23;
   this->rho = rho;
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticOrthotropic<3>>::testComputeStress() {
   UInt Dim = 3;
 
   // material frame of reference is rotate by rotation_matrix starting from
   // canonical basis
   Matrix<Real> rotation_matrix = getRandomRotation();
 
   // canonical basis as expressed in the material frame of reference, as
   // required by MaterialElasticOrthotropic class (it is simply given by the
   // columns of the rotation_matrix; the lines give the material basis expressed
   // in the canonical frame of reference)
   *this->dir_vecs[0] = rotation_matrix(0);
   *this->dir_vecs[1] = rotation_matrix(1);
   *this->dir_vecs[2] = rotation_matrix(2);
 
   // set internal Cijkl matrix expressed in the canonical frame of reference
   this->updateInternalParameters();
 
   // gradient in material frame of reference
   auto grad_u = this->getComposedStrain(2.);
 
   // gradient in canonical basis (we need to rotate *back* to the canonical
   // basis)
   auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);
 
   // stress in the canonical basis
   Matrix<Real> sigma_rot(3, 3);
   this->computeStressOnQuad(grad_u_rot, sigma_rot);
 
   // stress in the material reference (we need to apply the rotation)
   auto sigma = this->applyRotation(sigma_rot, rotation_matrix);
 
   // construction of Cijkl engineering tensor in the *material* frame of
   // reference
   // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
   Real nu21 = nu12 * E2 / E1;
   Real nu31 = nu13 * E3 / E1;
   Real nu32 = nu23 * E3 / E2;
   Real gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 -
                     2 * nu21 * nu32 * nu13);
 
   Matrix<Real> C_expected(6, 6);
   C_expected(0, 0) = gamma * E1 * (1 - nu23 * nu32);
   C_expected(1, 1) = gamma * E2 * (1 - nu13 * nu31);
   C_expected(2, 2) = gamma * E3 * (1 - nu12 * nu21);
 
   C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * (nu21 + nu31 * nu23);
   C_expected(2, 0) = C_expected(0, 2) = gamma * E1 * (nu31 + nu21 * nu32);
   C_expected(2, 1) = C_expected(1, 2) = gamma * E2 * (nu32 + nu12 * nu31);
 
   C_expected(3, 3) = G23;
   C_expected(4, 4) = G13;
   C_expected(5, 5) = G12;
 
   // epsilon is computed directly in the *material* frame of reference
   Matrix<Real> epsilon = 0.5 * (grad_u + grad_u.transpose());
 
   // sigma_expected is computed directly in the *material* frame of reference
   Matrix<Real> sigma_expected(Dim, Dim);
   for (UInt i = 0; i < Dim; ++i) {
     for (UInt j = 0; j < Dim; ++j) {
       sigma_expected(i, i) += C_expected(i, j) * epsilon(j, j);
     }
   }
 
   sigma_expected(0, 1) = C_expected(5, 5) * 2 * epsilon(0, 1);
   sigma_expected(0, 2) = C_expected(4, 4) * 2 * epsilon(0, 2);
   sigma_expected(1, 2) = C_expected(3, 3) * 2 * epsilon(1, 2);
   sigma_expected(1, 0) = sigma_expected(0, 1);
   sigma_expected(2, 0) = sigma_expected(0, 2);
   sigma_expected(2, 1) = sigma_expected(1, 2);
 
   // sigmas are checked in the *material* frame of reference
   auto diff = sigma - sigma_expected;
   Real stress_error = diff.norm<L_inf>();
   EXPECT_NEAR(stress_error, 0., 1e-13);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticOrthotropic<3>>::testEnergyDensity() {
   Matrix<Real> sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
   Matrix<Real> eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
   Real epot = 0;
   Real solution = 5.5;
   this->computePotentialEnergyOnQuad(eps, sigma, epot);
   EXPECT_NEAR(epot, solution, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticOrthotropic<3>>::testComputeTangentModuli() {
   // Note: for this test material and canonical basis coincide
   UInt Dim = 3;
 
   // construction of Cijkl engineering tensor in the *material* frame of
   // reference
   // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
   Real nu21 = nu12 * E2 / E1;
   Real nu31 = nu13 * E3 / E1;
   Real nu32 = nu23 * E3 / E2;
   Real gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 -
                     2 * nu21 * nu32 * nu13);
 
   Matrix<Real> C_expected(2 * Dim, 2 * Dim, 0);
   C_expected(0, 0) = gamma * E1 * (1 - nu23 * nu32);
   C_expected(1, 1) = gamma * E2 * (1 - nu13 * nu31);
   C_expected(2, 2) = gamma * E3 * (1 - nu12 * nu21);
 
   C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * (nu21 + nu31 * nu23);
   C_expected(2, 0) = C_expected(0, 2) = gamma * E1 * (nu31 + nu21 * nu32);
   C_expected(2, 1) = C_expected(1, 2) = gamma * E2 * (nu32 + nu12 * nu31);
 
   C_expected(3, 3) = G23;
   C_expected(4, 4) = G13;
   C_expected(5, 5) = G12;
 
   Matrix<Real> tangent(6, 6);
   this->computeTangentModuliOnQuad(tangent);
 
   Real tangent_error = (tangent - C_expected).norm<L_2>();
   EXPECT_NEAR(tangent_error, 0, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testCelerity() {
   // Note: for this test material and canonical basis coincide
   UInt Dim = 3;
   // construction of Cijkl engineering tensor in the *material* frame of
   // reference
   // ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
   Real nu21 = nu12 * E2 / E1;
   Real nu31 = nu13 * E3 / E1;
   Real nu32 = nu23 * E3 / E2;
   Real gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 -
                     2 * nu21 * nu32 * nu13);
 
   Matrix<Real> C_expected(2 * Dim, 2 * Dim, 0);
   C_expected(0, 0) = gamma * E1 * (1 - nu23 * nu32);
   C_expected(1, 1) = gamma * E2 * (1 - nu13 * nu31);
   C_expected(2, 2) = gamma * E3 * (1 - nu12 * nu21);
 
   C_expected(1, 0) = C_expected(0, 1) = gamma * E1 * (nu21 + nu31 * nu23);
   C_expected(2, 0) = C_expected(0, 2) = gamma * E1 * (nu31 + nu21 * nu32);
   C_expected(2, 1) = C_expected(1, 2) = gamma * E2 * (nu32 + nu12 * nu31);
 
   C_expected(3, 3) = G23;
   C_expected(4, 4) = G13;
   C_expected(5, 5) = G12;
 
   Vector<Real> eig_expected(6);
   C_expected.eig(eig_expected);
 
   auto celerity_expected = std::sqrt(eig_expected(0) / rho);
 
   auto celerity = this->getCelerity(Element());
 
   EXPECT_NEAR(celerity_expected, celerity, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::setParams() {
   Matrix<Real> C = {
       {1.0, 0.3, 0.4},
       {0.3, 2.0, 0.1},
       {0.4, 0.1, 1.5},
   };
 
   for (auto i = 0u; i < C.rows(); ++i)
     for (auto j = 0u; j < C.cols(); ++j)
       this->Cprime(i, j) = C(i, j);
 
   this->rho = 2.7;
 
   // material frame of reference is rotate by rotation_matrix starting from
   // canonical basis
   Matrix<Real> rotation_matrix = getRandomRotation();
 
   // canonical basis as expressed in the material frame of reference, as
   // required by MaterialElasticLinearAnisotropic class (it is simply given by
   // the columns of the rotation_matrix; the lines give the material basis
   // expressed in the canonical frame of reference)
   *this->dir_vecs[0] = rotation_matrix(0);
   *this->dir_vecs[1] = rotation_matrix(1);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testComputeStress() {
   Matrix<Real> C = {
       {1.0, 0.3, 0.4},
       {0.3, 2.0, 0.1},
       {0.4, 0.1, 1.5},
   };
 
   Matrix<Real> rotation_matrix(2, 2);
 
   rotation_matrix(0) = *this->dir_vecs[0];
   rotation_matrix(1) = *this->dir_vecs[1];
 
   // gradient in material frame of reference
   auto grad_u = this->getComposedStrain(1.).block(0, 0, 2, 2);
 
   // gradient in canonical basis (we need to rotate *back* to the canonical
   // basis)
   auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);
 
   // stress in the canonical basis
   Matrix<Real> sigma_rot(2, 2);
   this->computeStressOnQuad(grad_u_rot, sigma_rot);
 
   // stress in the material reference (we need to apply the rotation)
   auto sigma = this->applyRotation(sigma_rot, rotation_matrix);
 
   // epsilon is computed directly in the *material* frame of reference
   Matrix<Real> epsilon = 0.5 * (grad_u + grad_u.transpose());
   Vector<Real> epsilon_voigt(3);
   epsilon_voigt(0) = epsilon(0, 0);
   epsilon_voigt(1) = epsilon(1, 1);
   epsilon_voigt(2) = 2 * epsilon(0, 1);
 
   // sigma_expected is computed directly in the *material* frame of reference
   Vector<Real> sigma_voigt = C * epsilon_voigt;
   Matrix<Real> sigma_expected(2, 2);
   sigma_expected(0, 0) = sigma_voigt(0);
   sigma_expected(1, 1) = sigma_voigt(1);
   sigma_expected(0, 1) = sigma_expected(1, 0) = sigma_voigt(2);
 
   // sigmas are checked in the *material* frame of reference
   auto diff = sigma - sigma_expected;
   Real stress_error = diff.norm<L_inf>();
   EXPECT_NEAR(stress_error, 0., 1e-13);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testEnergyDensity() {
   Matrix<Real> sigma = {{1, 2}, {2, 4}};
   Matrix<Real> eps = {{1, 0}, {0, 1}};
   Real epot = 0;
   Real solution = 2.5;
   this->computePotentialEnergyOnQuad(eps, sigma, epot);
   EXPECT_NEAR(epot, solution, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<
     MaterialElasticLinearAnisotropic<2>>::testComputeTangentModuli() {
   Matrix<Real> tangent(3, 3);
   this->computeTangentModuliOnQuad(tangent);
 
   Real tangent_error = (tangent - C).norm<L_2>();
   EXPECT_NEAR(tangent_error, 0, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testCelerity() {
   Vector<Real> eig_expected(3);
   C.eig(eig_expected);
 
   auto celerity_expected = std::sqrt(eig_expected(0) / this->rho);
   auto celerity = this->getCelerity(Element());
 
   EXPECT_NEAR(celerity_expected, celerity, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::setParams() {
   // Note: for this test material and canonical basis coincide
   Matrix<Real> C = {
       {1.0, 0.3, 0.4, 0.3, 0.2, 0.1}, {0.3, 2.0, 0.1, 0.2, 0.3, 0.2},
       {0.4, 0.1, 1.5, 0.1, 0.4, 0.3}, {0.3, 0.2, 0.1, 2.4, 0.1, 0.4},
       {0.2, 0.3, 0.4, 0.1, 0.9, 0.1}, {0.1, 0.2, 0.3, 0.4, 0.1, 1.2},
   };
 
   for (auto i = 0u; i < C.rows(); ++i)
     for (auto j = 0u; j < C.cols(); ++j)
       this->Cprime(i, j) = C(i, j);
 
   this->rho = 2.9;
 
   // material frame of reference is rotate by rotation_matrix starting from
   // canonical basis
   Matrix<Real> rotation_matrix = getRandomRotation();
 
   // canonical basis as expressed in the material frame of reference, as
   // required by MaterialElasticLinearAnisotropic class (it is simply given by
   // the columns of the rotation_matrix; the lines give the material basis
   // expressed in the canonical frame of reference)
   *this->dir_vecs[0] = rotation_matrix(0);
   *this->dir_vecs[1] = rotation_matrix(1);
   *this->dir_vecs[2] = rotation_matrix(2);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testComputeStress() {
   Matrix<Real> C = {
       {1.0, 0.3, 0.4, 0.3, 0.2, 0.1}, {0.3, 2.0, 0.1, 0.2, 0.3, 0.2},
       {0.4, 0.1, 1.5, 0.1, 0.4, 0.3}, {0.3, 0.2, 0.1, 2.4, 0.1, 0.4},
       {0.2, 0.3, 0.4, 0.1, 0.9, 0.1}, {0.1, 0.2, 0.3, 0.4, 0.1, 1.2},
   };
 
   Matrix<Real> rotation_matrix(3, 3);
 
   rotation_matrix(0) = *this->dir_vecs[0];
   rotation_matrix(1) = *this->dir_vecs[1];
   rotation_matrix(2) = *this->dir_vecs[2];
 
   // gradient in material frame of reference
   auto grad_u = this->getComposedStrain(2.);
 
   // gradient in canonical basis (we need to rotate *back* to the canonical
   // basis)
   auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);
 
   // stress in the canonical basis
   Matrix<Real> sigma_rot(3, 3);
   this->computeStressOnQuad(grad_u_rot, sigma_rot);
 
   // stress in the material reference (we need to apply the rotation)
   auto sigma = this->applyRotation(sigma_rot, rotation_matrix);
 
   // epsilon is computed directly in the *material* frame of reference
   Matrix<Real> epsilon = 0.5 * (grad_u + grad_u.transpose());
   Vector<Real> epsilon_voigt(6);
   epsilon_voigt(0) = epsilon(0, 0);
   epsilon_voigt(1) = epsilon(1, 1);
   epsilon_voigt(2) = epsilon(2, 2);
   epsilon_voigt(3) = 2 * epsilon(1, 2);
   epsilon_voigt(4) = 2 * epsilon(0, 2);
   epsilon_voigt(5) = 2 * epsilon(0, 1);
 
   // sigma_expected is computed directly in the *material* frame of reference
   Vector<Real> sigma_voigt = C * epsilon_voigt;
   Matrix<Real> sigma_expected(3, 3);
   sigma_expected(0, 0) = sigma_voigt(0);
   sigma_expected(1, 1) = sigma_voigt(1);
   sigma_expected(2, 2) = sigma_voigt(2);
   sigma_expected(1, 2) = sigma_expected(2, 1) = sigma_voigt(3);
   sigma_expected(0, 2) = sigma_expected(2, 0) = sigma_voigt(4);
   sigma_expected(0, 1) = sigma_expected(1, 0) = sigma_voigt(5);
 
   // sigmas are checked in the *material* frame of reference
   auto diff = sigma - sigma_expected;
   Real stress_error = diff.norm<L_inf>();
   EXPECT_NEAR(stress_error, 0., 1e-13);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testEnergyDensity() {
   Matrix<Real> sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
   Matrix<Real> eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
   Real epot = 0;
   Real solution = 5.5;
   this->computePotentialEnergyOnQuad(eps, sigma, epot);
   EXPECT_NEAR(epot, solution, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<
     MaterialElasticLinearAnisotropic<3>>::testComputeTangentModuli() {
   Matrix<Real> tangent(6, 6);
   this->computeTangentModuliOnQuad(tangent);
 
   Real tangent_error = (tangent - C).norm<L_2>();
   EXPECT_NEAR(tangent_error, 0, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testCelerity() {
   Vector<Real> eig_expected(6);
   C.eig(eig_expected);
 
   auto celerity_expected = std::sqrt(eig_expected(0) / this->rho);
 
   auto celerity = this->getCelerity(Element());
 
   EXPECT_NEAR(celerity_expected, celerity, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 
 namespace {
 
 template <typename T>
 class TestElasticMaterialFixture : public ::TestMaterialFixture<T> {};
 
-TYPED_TEST_SUITE(TestElasticMaterialFixture, mat_types);
+TYPED_TEST_SUITE(TestElasticMaterialFixture, mat_types, );
 
 TYPED_TEST(TestElasticMaterialFixture, ComputeStress) {
   this->material->testComputeStress();
 }
 
 TYPED_TEST(TestElasticMaterialFixture, EnergyDensity) {
   this->material->testEnergyDensity();
 }
 
 TYPED_TEST(TestElasticMaterialFixture, ComputeTangentModuli) {
   this->material->testComputeTangentModuli();
 }
 
 TYPED_TEST(TestElasticMaterialFixture, ComputeCelerity) {
   this->material->testCelerity();
 }
 
 } // namespace
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_finite_def_materials.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_finite_def_materials.cc
index 917026c44..a93c6a31c 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_finite_def_materials.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_finite_def_materials.cc
@@ -1,86 +1,86 @@
 /**
  * @file   test_finite_def_materials.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @date creation: Fri Nov 17 2017
  * @date last modification: Tue Feb 20 2018
  *
  * @brief
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "test_gtest_utils.hh"
 #include "test_material_fixtures.hh"
 /* -------------------------------------------------------------------------- */
 #include <material_neohookean.hh>
 #include <solid_mechanics_model.hh>
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 using mat_types = ::testing::Types<Traits<MaterialNeohookean, 1>,
                                    Traits<MaterialNeohookean, 2>,
                                    Traits<MaterialNeohookean, 3>>;
 
 /*****************************************************************/
 
 template <> void FriendMaterial<MaterialNeohookean<3>>::testComputeStress() {
   AKANTU_TO_IMPLEMENT();
 }
 
 /*****************************************************************/
 template <>
 void FriendMaterial<MaterialNeohookean<3>>::testComputeTangentModuli() {
   AKANTU_TO_IMPLEMENT();
 }
 
 /*****************************************************************/
 
 template <> void FriendMaterial<MaterialNeohookean<3>>::testEnergyDensity() {
   AKANTU_TO_IMPLEMENT();
 }
 
 /*****************************************************************/
 
 namespace {
 
 template <typename T>
 class TestFiniteDefMaterialFixture : public ::TestMaterialFixture<T> {};
 
-TYPED_TEST_SUITE(TestFiniteDefMaterialFixture, mat_types);
+TYPED_TEST_SUITE(TestFiniteDefMaterialFixture, mat_types, );
 
 TYPED_TEST(TestFiniteDefMaterialFixture, DISABLED_ComputeStress) {
   this->material->testComputeStress();
 }
 TYPED_TEST(TestFiniteDefMaterialFixture, DISABLED_EnergyDensity) {
   this->material->testEnergyDensity();
 }
 TYPED_TEST(TestFiniteDefMaterialFixture, DISABLED_ComputeTangentModuli) {
   this->material->testComputeTangentModuli();
 }
 TYPED_TEST(TestFiniteDefMaterialFixture, DISABLED_DefComputeCelerity) {
   this->material->testCelerity();
 }
 } // namespace
 /*****************************************************************/
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_thermal.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_thermal.cc
index 648ee7388..c84d43d8a 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_thermal.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_thermal.cc
@@ -1,104 +1,104 @@
 /**
  * @file   test_material_thermal.cc
  *
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  *
  * @date creation: Wed Aug 04 2010
  * @date last modification: Mon Jan 29 2018
  *
  * @brief  Test the thermal material
  *
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "material_thermal.hh"
 #include "solid_mechanics_model.hh"
 #include "test_material_fixtures.hh"
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 using mat_types =
     ::testing::Types<Traits<MaterialThermal, 1>, Traits<MaterialThermal, 2>,
                      Traits<MaterialThermal, 3>>;
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialThermal<3>>::testComputeStress() {
   Real E = 1.;
   Real nu = .3;
   Real alpha = 2;
   setParam("E", E);
   setParam("nu", nu);
   setParam("alpha", alpha);
 
   Real deltaT = 1;
   Real sigma = 0;
   this->computeStressOnQuad(sigma, deltaT);
   Real solution = -E / (1 - 2 * nu) * alpha * deltaT;
   auto error = std::abs(sigma - solution);
   ASSERT_NEAR(error, 0, 1e-14);
 }
 
 template <> void FriendMaterial<MaterialThermal<2>>::testComputeStress() {
   Real E = 1.;
   Real nu = .3;
   Real alpha = 2;
   setParam("E", E);
   setParam("nu", nu);
   setParam("alpha", alpha);
 
   Real deltaT = 1;
   Real sigma = 0;
   this->computeStressOnQuad(sigma, deltaT);
   Real solution = -E / (1 - 2 * nu) * alpha * deltaT;
   auto error = std::abs(sigma - solution);
   ASSERT_NEAR(error, 0, 1e-14);
 }
 
 template <> void FriendMaterial<MaterialThermal<1>>::testComputeStress() {
   Real E = 1.;
   Real nu = .3;
   Real alpha = 2;
   setParam("E", E);
   setParam("nu", nu);
   setParam("alpha", alpha);
 
   Real deltaT = 1;
   Real sigma = 0;
   this->computeStressOnQuad(sigma, deltaT);
   Real solution = -E * alpha * deltaT;
   auto error = std::abs(sigma - solution);
   ASSERT_NEAR(error, 0, 1e-14);
 }
 
 namespace {
 
 template <typename T>
 class TestMaterialThermalFixture : public ::TestMaterialFixture<T> {};
 
-TYPED_TEST_SUITE(TestMaterialThermalFixture, mat_types);
+TYPED_TEST_SUITE(TestMaterialThermalFixture, mat_types, );
 
 TYPED_TEST(TestMaterialThermalFixture, ThermalComputeStress) {
   this->material->testComputeStress();
 }
 } // namespace
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_plastic_materials.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_plastic_materials.cc
index dec293699..dd1ed42a9 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_plastic_materials.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_plastic_materials.cc
@@ -1,191 +1,191 @@
 /**
  * @file   test_plastic_materials.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @date creation: Fri Nov 17 2017
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  Tests the plastic material
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "material_linear_isotropic_hardening.hh"
 #include "solid_mechanics_model.hh"
 #include "test_material_fixtures.hh"
 #include <gtest/gtest.h>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 using mat_types = ::testing::Types<
     // Traits<MaterialLinearIsotropicHardening, 1>,
     // Traits<MaterialLinearIsotropicHardening, 2>,
     Traits<MaterialLinearIsotropicHardening, 3>>;
 
 /* -------------------------------------------------------------------------- */
 
 template <>
 void FriendMaterial<MaterialLinearIsotropicHardening<3>>::testComputeStress() {
 
   Real E = 1.;
   // Real nu = .3;
   Real nu = 0.;
   Real rho = 1.;
   Real sigma_0 = 1.;
   Real h = 0.;
   Real bulk_modulus_K = E / 3. / (1 - 2. * nu);
   Real shear_modulus_mu = 0.5 * E / (1 + nu);
 
   setParam("E", E);
   setParam("nu", nu);
   setParam("rho", rho);
   setParam("sigma_y", sigma_0);
   setParam("h", h);
 
   auto rotation_matrix = getRandomRotation();
 
   Real max_strain = 10.;
   Real strain_steps = 100;
   Real dt = max_strain / strain_steps;
   std::vector<double> steps(strain_steps);
   std::iota(steps.begin(), steps.end(), 0.);
 
   Matrix<Real> previous_grad_u_rot(3, 3, 0.);
   Matrix<Real> previous_sigma(3, 3, 0.);
   Matrix<Real> previous_sigma_rot(3, 3, 0.);
   Matrix<Real> inelastic_strain_rot(3, 3, 0.);
   Matrix<Real> inelastic_strain(3, 3, 0.);
   Matrix<Real> previous_inelastic_strain(3, 3, 0.);
   Matrix<Real> previous_inelastic_strain_rot(3, 3, 0.);
   Matrix<Real> sigma_rot(3, 3, 0.);
   Real iso_hardening = 0.;
   Real previous_iso_hardening = 0.;
 
   // hydrostatic loading (should not plastify)
   for (auto && i : steps) {
     auto t = i * dt;
 
     auto grad_u = this->getHydrostaticStrain(t);
     auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix);
 
     this->computeStressOnQuad(grad_u_rot, previous_grad_u_rot, sigma_rot,
                               previous_sigma_rot, inelastic_strain_rot,
                               previous_inelastic_strain_rot, iso_hardening,
                               previous_iso_hardening, 0., 0.);
 
     auto sigma = this->reverseRotation(sigma_rot, rotation_matrix);
 
     Matrix<Real> sigma_expected =
         t * 3. * bulk_modulus_K * Matrix<Real>::eye(3, 1.);
 
     Real stress_error = (sigma - sigma_expected).norm<L_inf>();
 
     ASSERT_NEAR(stress_error, 0., 1e-13);
     ASSERT_NEAR(inelastic_strain_rot.norm<L_inf>(), 0., 1e-13);
 
     previous_grad_u_rot = grad_u_rot;
     previous_sigma_rot = sigma_rot;
     previous_inelastic_strain_rot = inelastic_strain_rot;
     previous_iso_hardening = iso_hardening;
   }
 
   // deviatoric loading (should plastify)
   // stress at onset of plastication
   Real beta = sqrt(42);
   Real t_P = sigma_0 / 2. / shear_modulus_mu / beta;
   Matrix<Real> sigma_P = sigma_0 / beta * this->getDeviatoricStrain(1.);
 
   for (auto && i : steps) {
 
     auto t = i * dt;
     auto grad_u = this->getDeviatoricStrain(t);
     auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix);
     Real iso_hardening, previous_iso_hardening;
 
     this->computeStressOnQuad(grad_u_rot, previous_grad_u_rot, sigma_rot,
                               previous_sigma_rot, inelastic_strain_rot,
                               previous_inelastic_strain_rot, iso_hardening,
                               previous_iso_hardening, 0., 0.);
 
     auto sigma = this->reverseRotation(sigma_rot, rotation_matrix);
     auto inelastic_strain =
         this->reverseRotation(inelastic_strain_rot, rotation_matrix);
 
     if (t < t_P) {
 
       Matrix<Real> sigma_expected =
           shear_modulus_mu * (grad_u + grad_u.transpose());
 
       Real stress_error = (sigma - sigma_expected).norm<L_inf>();
       ASSERT_NEAR(stress_error, 0., 1e-13);
       ASSERT_NEAR(inelastic_strain_rot.norm<L_inf>(), 0., 1e-13);
     } else if (t > t_P + dt) {
       // skip the transition from non plastic to plastic
 
       auto delta_lambda_expected =
           dt / t * previous_sigma.doubleDot(grad_u + grad_u.transpose()) / 2.;
       auto delta_inelastic_strain_expected =
           delta_lambda_expected * 3. / 2. / sigma_0 * previous_sigma;
       auto inelastic_strain_expected =
           delta_inelastic_strain_expected + previous_inelastic_strain;
       ASSERT_NEAR((inelastic_strain - inelastic_strain_expected).norm<L_inf>(),
                   0., 1e-13);
       auto delta_sigma_expected =
           2. * shear_modulus_mu *
           (0.5 * dt / t * (grad_u + grad_u.transpose()) -
            delta_inelastic_strain_expected);
 
       auto delta_sigma = sigma - previous_sigma;
       ASSERT_NEAR((delta_sigma_expected - delta_sigma).norm<L_inf>(), 0.,
                   1e-13);
     }
     previous_sigma = sigma;
     previous_sigma_rot = sigma_rot;
     previous_grad_u_rot = grad_u_rot;
     previous_inelastic_strain = inelastic_strain;
     previous_inelastic_strain_rot = inelastic_strain_rot;
   }
 }
 
 namespace {
 
 template <typename T>
 class TestPlasticMaterialFixture : public ::TestMaterialFixture<T> {};
 
-TYPED_TEST_SUITE(TestPlasticMaterialFixture, mat_types);
+TYPED_TEST_SUITE(TestPlasticMaterialFixture, mat_types, );
 
 TYPED_TEST(TestPlasticMaterialFixture, ComputeStress) {
   this->material->testComputeStress();
 }
 TYPED_TEST(TestPlasticMaterialFixture, DISABLED_EnergyDensity) {
   this->material->testEnergyDensity();
 }
 TYPED_TEST(TestPlasticMaterialFixture, DISABLED_ComputeTangentModuli) {
   this->material->testComputeTangentModuli();
 }
 TYPED_TEST(TestPlasticMaterialFixture, DISABLED_ComputeCelerity) {
   this->material->testCelerity();
 }
 } // namespace
 
 /*****************************************************************/
diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc
index 61fdad515..441dbdc00 100644
--- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc
+++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc
@@ -1,320 +1,320 @@
 /**
  * @file   test_solid_mechanics_model_dynamics.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @date creation: Wed Nov 29 2017
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  test of the class SolidMechanicsModel on the 3d cube
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "boundary_condition_functor.hh"
 #include "test_solid_mechanics_model_fixture.hh"
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 namespace {
 
 const Real A = 1e-1;
 const Real E = 1.;
 const Real poisson = 3. / 10;
 const Real lambda = E * poisson / (1 + poisson) / (1 - 2 * poisson);
 const Real mu = E / 2 / (1. + poisson);
 const Real rho = 1.;
 const Real cp = std::sqrt((lambda + 2 * mu) / rho);
 const Real cs = std::sqrt(mu / rho);
 const Real c = std::sqrt(E / rho);
 
 const Vector<Real> k = {.5, 0., 0.};
 const Vector<Real> psi1 = {0., 0., 1.};
 const Vector<Real> psi2 = {0., 1., 0.};
 const Real knorm = k.norm();
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim> struct Verification {};
 
 /* -------------------------------------------------------------------------- */
 template <> struct Verification<1> {
   void displacement(Vector<Real> & disp, const Vector<Real> & coord,
                     Real current_time) {
     const auto & x = coord(_x);
     const Real omega = c * k[0];
     disp(_x) = A * cos(k[0] * x - omega * current_time);
   }
 
   void velocity(Vector<Real> & vel, const Vector<Real> & coord,
                 Real current_time) {
     const auto & x = coord(_x);
     const Real omega = c * k[0];
     vel(_x) = omega * A * sin(k[0] * x - omega * current_time);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <> struct Verification<2> {
   void displacement(Vector<Real> & disp, const Vector<Real> & X,
                     Real current_time) {
     Vector<Real> kshear = {k[1], k[0]};
     Vector<Real> kpush = {k[0], k[1]};
 
     const Real omega_p = knorm * cp;
     const Real omega_s = knorm * cs;
 
     Real phase_p = X.dot(kpush) - omega_p * current_time;
     Real phase_s = X.dot(kpush) - omega_s * current_time;
 
     disp = A * (kpush * cos(phase_p) + kshear * cos(phase_s));
   }
 
   void velocity(Vector<Real> & vel, const Vector<Real> & X, Real current_time) {
     Vector<Real> kshear = {k[1], k[0]};
     Vector<Real> kpush = {k[0], k[1]};
 
     const Real omega_p = knorm * cp;
     const Real omega_s = knorm * cs;
 
     Real phase_p = X.dot(kpush) - omega_p * current_time;
     Real phase_s = X.dot(kpush) - omega_s * current_time;
 
     vel =
         A * (kpush * omega_p * cos(phase_p) + kshear * omega_s * cos(phase_s));
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <> struct Verification<3> {
   void displacement(Vector<Real> & disp, const Vector<Real> & coord,
                     Real current_time) {
     const auto & X = coord;
     Vector<Real> kpush = k;
     Vector<Real> kshear1(3);
     Vector<Real> kshear2(3);
     kshear1.crossProduct(k, psi1);
     kshear2.crossProduct(k, psi2);
 
     const Real omega_p = knorm * cp;
     const Real omega_s = knorm * cs;
 
     Real phase_p = X.dot(kpush) - omega_p * current_time;
     Real phase_s = X.dot(kpush) - omega_s * current_time;
 
     disp = A * (kpush * cos(phase_p) + kshear1 * cos(phase_s) +
                 kshear2 * cos(phase_s));
   }
 
   void velocity(Vector<Real> & vel, const Vector<Real> & coord,
                 Real current_time) {
     const auto & X = coord;
     Vector<Real> kpush = k;
     Vector<Real> kshear1(3);
     Vector<Real> kshear2(3);
     kshear1.crossProduct(k, psi1);
     kshear2.crossProduct(k, psi2);
 
     const Real omega_p = knorm * cp;
     const Real omega_s = knorm * cs;
 
     Real phase_p = X.dot(kpush) - omega_p * current_time;
     Real phase_s = X.dot(kpush) - omega_s * current_time;
 
     vel =
         A * (kpush * omega_p * cos(phase_p) + kshear1 * omega_s * cos(phase_s) +
              kshear2 * omega_s * cos(phase_s));
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <ElementType _type>
 class SolutionFunctor : public BC::Dirichlet::DirichletFunctor {
 public:
   SolutionFunctor(Real current_time, SolidMechanicsModel & model)
       : current_time(current_time), model(model) {}
 
 public:
   static constexpr UInt dim = ElementClass<_type>::getSpatialDimension();
 
   inline void operator()(UInt node, Vector<bool> & flags, Vector<Real> & primal,
                          const Vector<Real> & coord) const {
     flags(0) = true;
     auto & vel = model.getVelocity();
     auto it = vel.begin(model.getSpatialDimension());
     Vector<Real> v = it[node];
 
     Verification<dim> verif;
     verif.displacement(primal, coord, current_time);
     verif.velocity(v, coord, current_time);
   }
 
 private:
   Real current_time;
   SolidMechanicsModel & model;
 };
 
 /* -------------------------------------------------------------------------- */
 // This fixture uses somewhat finer meshes representing bars.
 template <typename type_, typename analysis_method_>
 class TestSMMFixtureBar : public TestSMMFixture<type_> {
   using parent = TestSMMFixture<type_>;
 
 public:
   void SetUp() override {
     this->mesh_file =
         "../patch_tests/data/bar" + std::to_string(this->type) + ".msh";
     parent::SetUp();
 
     auto analysis_method = analysis_method_::value;
     this->initModel("test_solid_mechanics_model_"
                     "dynamics_material.dat",
                     analysis_method);
 
     const auto & position = this->mesh->getNodes();
     auto & displacement = this->model->getDisplacement();
     auto & velocity = this->model->getVelocity();
 
     constexpr auto dim = parent::spatial_dimension;
 
     Verification<dim> verif;
 
     for (auto && tuple :
          zip(make_view(position, dim), make_view(displacement, dim),
              make_view(velocity, dim))) {
       verif.displacement(std::get<1>(tuple), std::get<0>(tuple), 0.);
       verif.velocity(std::get<2>(tuple), std::get<0>(tuple), 0.);
     }
 
     if (this->dump_paraview)
       this->model->dump();
 
     /// boundary conditions
     this->model->applyBC(SolutionFunctor<parent::type>(0., *this->model), "BC");
 
     wave_velocity = 1.; // sqrt(E/rho) = sqrt(1/1) = 1
     simulation_time = 5 / wave_velocity;
     time_step = this->model->getTimeStep();
 
     max_steps = simulation_time / time_step; // 100
   }
 
   void solveStep() {
     constexpr auto dim = parent::spatial_dimension;
     Real current_time = 0.;
     const auto & position = this->mesh->getNodes();
     const auto & displacement = this->model->getDisplacement();
     UInt nb_nodes = this->mesh->getNbNodes();
     UInt nb_global_nodes = this->mesh->getNbGlobalNodes();
 
     max_error = -1.;
 
     Array<Real> displacement_solution(nb_nodes, dim);
     Verification<dim> verif;
 
     auto ndump = 50;
     auto dump_freq = max_steps / ndump;
 
     for (UInt s = 0; s < this->max_steps;
          ++s, current_time += this->time_step) {
       if (s % dump_freq == 0 && this->dump_paraview)
         this->model->dump();
 
       /// boundary conditions
       this->model->applyBC(
           SolutionFunctor<parent::type>(current_time, *this->model), "BC");
 
       // compute the disp solution
       for (auto && tuple : zip(make_view(position, dim),
                                make_view(displacement_solution, dim))) {
         verif.displacement(std::get<1>(tuple), std::get<0>(tuple),
                            current_time);
       }
 
       // compute the error solution
       Real disp_error = 0.;
       auto n = 0;
       for (auto && tuple : zip(make_view(displacement, dim),
                                make_view(displacement_solution, dim))) {
         if (this->mesh->isLocalOrMasterNode(n)) {
           auto diff = std::get<1>(tuple) - std::get<0>(tuple);
           disp_error += diff.dot(diff);
         }
         ++n;
       }
 
       this->mesh->getCommunicator().allReduce(disp_error,
                                               SynchronizerOperation::_sum);
 
       disp_error = sqrt(disp_error) / nb_global_nodes;
       max_error = std::max(disp_error, max_error);
 
       this->model->solveStep();
     }
   }
 
 protected:
   Real time_step;
   Real wave_velocity;
   Real simulation_time;
   UInt max_steps;
   Real max_error{-1};
 };
 
 template <AnalysisMethod t>
 using analysis_method_t = std::integral_constant<AnalysisMethod, t>;
 
 using TestTypes = gtest_list_t<TestElementTypes>;
 
 template <typename type_>
 using TestSMMFixtureBarExplicit =
     TestSMMFixtureBar<type_, analysis_method_t<_explicit_lumped_mass>>;
 
-TYPED_TEST_SUITE(TestSMMFixtureBarExplicit, TestTypes);
+TYPED_TEST_SUITE(TestSMMFixtureBarExplicit, TestTypes, );
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(TestSMMFixtureBarExplicit, Dynamics) {
   this->solveStep();
   EXPECT_NEAR(this->max_error, 0., 2e-3);
   // std::cout << "max error: " << max_error << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_IMPLICIT)
 template <typename type_>
 using TestSMMFixtureBarImplicit =
     TestSMMFixtureBar<type_, analysis_method_t<_implicit_dynamic>>;
 
-TYPED_TEST_SUITE(TestSMMFixtureBarImplicit, TestTypes);
+TYPED_TEST_SUITE(TestSMMFixtureBarImplicit, TestTypes, );
 
 TYPED_TEST(TestSMMFixtureBarImplicit, Dynamics) {
   if (this->type == _segment_2 and
       (this->mesh->getCommunicator().getNbProc() > 2)) {
     // The error are just to high after (hopefully because of the two small test
     // case)
     SUCCEED();
     return;
   }
   this->solveStep();
   EXPECT_NEAR(this->max_error, 0., 2e-3);
 }
 #endif
 
 } // namespace
diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh
index 3f2cd46fb..3b6e5e232 100644
--- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh
+++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh
@@ -1,129 +1,129 @@
 /**
  * @file   test_solid_mechanics_model_fixture.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Nov 14 2017
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Main solif mechanics test file
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "communicator.hh"
 #include "solid_mechanics_model.hh"
 #include "test_gtest_utils.hh"
 #include "mesh_utils.hh"
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_TEST_SOLID_MECHANICS_MODEL_FIXTURE_HH_
 #define AKANTU_TEST_SOLID_MECHANICS_MODEL_FIXTURE_HH_
 
 using namespace akantu;
 
 // This fixture uses very small meshes with a volume of 1.
 template <typename type_> class TestSMMFixture : public ::testing::Test {
 public:
   static constexpr ElementType type = type_::value;
   static constexpr size_t spatial_dimension =
       ElementClass<type>::getSpatialDimension();
 
   void SetUp() override {
     this->mesh = std::make_unique<Mesh>(this->spatial_dimension);
     auto prank = Communicator::getStaticCommunicator().whoAmI();
     if (prank == 0) {
       this->mesh->read(this->mesh_file);
       if(spatial_dimension > 1 and mesh->getNbElement(spatial_dimension - 1) == 0) {
         MeshUtils::buildFacets(*this->mesh);
       }
     }
 
     mesh->distribute();
 
     SCOPED_TRACE(std::to_string(this->type).c_str());
 
     model = std::make_unique<SolidMechanicsModel>(*mesh, _all_dimensions,
                                                   std::to_string(this->type));
   }
 
   void initModel(const ID & input, const AnalysisMethod & analysis_method) {
     getStaticParser().parse(input);
     this->model->initFull(_analysis_method = analysis_method);
 
     if (analysis_method != _static) {
       auto time_step = this->model->getStableTimeStep() / 10.;
       this->model->setTimeStep(time_step);
     }
     // std::cout << "timestep: " << time_step << std::endl;
 
     if (this->dump_paraview) {
       std::stringstream base_name;
       base_name << "bar" << analysis_method << this->type;
       this->model->setBaseName(base_name.str());
       this->model->addDumpFieldVector("displacement");
       this->model->addDumpFieldVector("blocked_dofs");
       if (analysis_method != _static) {
         this->model->addDumpField("velocity");
         this->model->addDumpField("acceleration");
       }
       if (this->mesh->isDistributed()) {
         this->model->addDumpField("partitions");
       }
       this->model->addDumpFieldVector("external_force");
       this->model->addDumpFieldVector("internal_force");
       this->model->addDumpField("stress");
       this->model->addDumpField("strain");
     }
   }
 
   void TearDown() override {
     model.reset(nullptr);
     mesh.reset(nullptr);
   }
 
 protected:
   std::string mesh_file{std::to_string(this->type) + ".msh"};
   std::unique_ptr<Mesh> mesh;
   std::unique_ptr<SolidMechanicsModel> model;
   bool dump_paraview{false};
 };
 
 template <typename type_> constexpr ElementType TestSMMFixture<type_>::type;
 template <typename type_>
 constexpr size_t TestSMMFixture<type_>::spatial_dimension;
 
 template <typename T>
 using is_not_pentahedron =
     aka::negation<aka::disjunction<is_element<T, _pentahedron_6>,
                                    is_element<T, _pentahedron_15>>>;
 
 using TestElementTypesFiltered =
     tuple_filter_t<is_not_pentahedron, TestElementTypes>;
 
 // using gtest_element_types = gtest_list_t<TestElementTypesFiltered>;
 using gtest_element_types = gtest_list_t<TestElementTypes>;
 
-TYPED_TEST_SUITE(TestSMMFixture, gtest_element_types);
+TYPED_TEST_SUITE(TestSMMFixture, gtest_element_types, );
 
 #endif /* AKANTU_TEST_SOLID_MECHANICS_MODEL_FIXTURE_HH_ */