diff --git a/packages/eigen.cmake b/packages/eigen.cmake
index 4733a342a..11f4c2e80 100644
--- a/packages/eigen.cmake
+++ b/packages/eigen.cmake
@@ -1,12 +1,12 @@
 package_declare(Eigen3 EXTERNAL
   DESCRIPTION "Add Eigen3 dependency to akantu"
   SYSTEM AUTO third-party/cmake/eigen3.cmake
-  EXTRA_PACKAGE_OPTIONS ARGS 3.3
+  EXTRA_PACKAGE_OPTIONS ARGS 3.4
   COMPILE_FLAGS CXX -DEIGEN_MAX_CPP_VER=14
   )
 
 mark_as_advanced(Eigen3_DIR)
 package_add_third_party_script_variable(Eigen3
-  EIGEN3_VERSION "3.3.9")
+  EIGEN3_VERSION "3.4.0")
 package_add_third_party_script_variable(Eigen3
   EIGEN3_GIT "https://gitlab.com/libeigen/eigen.git")
diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
index 4c37449bf..b70afd833 100644
--- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
+++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
@@ -1,269 +1,269 @@
 /**
  * @file   material_elastic_linear_anisotropic.cc
  *
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Till Junge <till.junge@epfl.ch>
  * @author Enrico Milanese <enrico.milanese@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Sep 25 2013
  * @date last modification: Fri Jul 24 2020
  *
  * @brief  Anisotropic elastic material
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2014-2021 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_elastic_linear_anisotropic.hh"
 #include "solid_mechanics_model.hh"
 #include <algorithm>
 #include <sstream>
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 template <Int dim>
 MaterialElasticLinearAnisotropic<dim>::MaterialElasticLinearAnisotropic(
     SolidMechanicsModel & model, const ID & id, bool symmetric)
     : Material(model, id), rot_mat(dim, dim), Cprime(dim * dim, dim * dim),
       C(voigt_h::size, voigt_h::size), eigC(voigt_h::size),
       symmetric(symmetric), was_stiffness_assembled(false) {
   AKANTU_DEBUG_IN();
 
   this->dir_vecs.push_back(std::make_unique<Vector<Real, dim>>());
   (*this->dir_vecs.back())[0] = 1.;
   this->registerParam("n1", *(this->dir_vecs.back()), _pat_parsmod,
                       "Direction of main material axis");
 
   if (dim > 1) {
     this->dir_vecs.push_back(std::make_unique<Vector<Real, dim>>());
     (*this->dir_vecs.back())[1] = 1.;
     this->registerParam("n2", *(this->dir_vecs.back()), _pat_parsmod,
                         "Direction of secondary material axis");
   }
 
   if (dim > 2) {
     this->dir_vecs.push_back(std::make_unique<Vector<Real, dim>>());
     (*this->dir_vecs.back())[2] = 1.;
     this->registerParam("n3", *(this->dir_vecs.back()), _pat_parsmod,
                         "Direction of tertiary material axis");
   }
 
   for (auto i : arange(voigt_h::size)) {
     decltype(i) start = 0;
     if (this->symmetric) {
       start = i;
     }
     for (auto j : arange(start, voigt_h::size)) {
       auto param = "C" + std::to_string(i + 1) + std::to_string(j + 1);
       this->registerParam(param, this->Cprime(i, j), Real(0.), _pat_parsmod,
                           "Coefficient " + param);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <Int dim> void MaterialElasticLinearAnisotropic<dim>::initMaterial() {
   AKANTU_DEBUG_IN();
   Material::initMaterial();
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <Int dim>
 void MaterialElasticLinearAnisotropic<dim>::updateInternalParameters() {
 
   Material::updateInternalParameters();
   if (this->symmetric) {
     for (auto i : arange(voigt_h::size)) {
       for (auto j : arange(i + 1, voigt_h::size)) {
         this->Cprime(j, i) = this->Cprime(i, j);
       }
     }
   }
   this->rotateCprime();
   this->C.eig(this->eigC);
 
   this->was_stiffness_assembled = false;
 }
 
 /* -------------------------------------------------------------------------- */
 template <Int Dim> void MaterialElasticLinearAnisotropic<Dim>::rotateCprime() {
   // start by filling the empty parts fo Cprime
   auto diff = Dim * Dim - voigt_h::size;
   for (auto i : arange(voigt_h::size, Dim * Dim)) {
     for (auto j : arange(Dim * Dim)) {
       this->Cprime(i, j) = this->Cprime(i - diff, j);
     }
   }
   for (auto i : arange(Dim * Dim)) {
     for (auto j : arange(voigt_h::size, Dim * Dim)) {
       this->Cprime(i, j) = this->Cprime(i, j - diff);
     }
   }
   // construction of rotator tensor
   // normalise rotation matrix
   for (auto j : arange(Dim)) {
     auto && rot_vec = this->rot_mat(j);
     rot_vec = *this->dir_vecs[j];
     rot_vec.normalize();
   }
 
   // make sure the vectors form a right-handed base
   Vector<Real, 3> v1, v2, v3;
   v3.zero();
 
   if (Dim == 2) {
     for (auto i : arange(Dim)) {
       v1[i] = this->rot_mat(0, i);
       v2[i] = this->rot_mat(1, i);
     }
 
     v3 = v1.cross(v2);
     if (v3.norm() < 8 * std::numeric_limits<Real>::epsilon()) {
       AKANTU_ERROR("The axis vectors parallel.");
     }
 
     v3.normalize();
   } else if (Dim == 3) {
     v1 = this->rot_mat(0);
     v2 = this->rot_mat(1);
     v3 = this->rot_mat(2);
   }
 
-  auto test_axis = v1.cross(v2) - v3;
+  Vector<Real, 3> test_axis = v1.cross(v2) - v3;
 
   if (test_axis.norm() > 8 * std::numeric_limits<Real>::epsilon()) {
     AKANTU_ERROR("The axis vectors do not form a right-handed coordinate "
                  << "system. I. e., ||n1 x n2 - n3|| should be zero, but "
                  << "it is " << test_axis.norm() << ".");
   }
 
   // create the rotator and the reverse rotator
   Matrix<Real, Dim * Dim, Dim * Dim> rotator;
   Matrix<Real, Dim * Dim, Dim * Dim> revrotor;
   for (auto i : arange(Dim)) {
     for (auto j : arange(Dim)) {
       for (auto k : arange(Dim)) {
         for (auto l : arange(Dim)) {
           auto I = voigt_h::mat[i][j];
           auto J = voigt_h::mat[k][l];
           rotator(I, J) = this->rot_mat(k, i) * this->rot_mat(l, j);
           revrotor(I, J) = this->rot_mat(i, k) * this->rot_mat(j, l);
         }
       }
     }
   }
 
   // create the full rotated matrix
   Matrix<Real, Dim * Dim, Dim * Dim> Cfull(Dim * Dim, Dim * Dim);
   Cfull = rotator * Cprime * revrotor;
 
   for (auto i : arange(voigt_h::size)) {
     for (auto j : arange(voigt_h::size)) {
       this->C(i, j) = Cfull(i, j);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <Int dim>
 void MaterialElasticLinearAnisotropic<dim>::computeStress(
     ElementType el_type, GhostType ghost_type) {
   auto && arguments = getArguments(el_type, ghost_type);
   for (auto && data : arguments) {
     this->computeStressOnQuad(data);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <Int dim>
 void MaterialElasticLinearAnisotropic<dim>::computeTangentModuli(
     ElementType el_type, Array<Real> & tangent_matrix, GhostType ghost_type) {
   auto && arguments =
       Material::getArgumentsTangent<dim>(tangent_matrix, el_type, ghost_type);
 
   for (auto && args : arguments) {
     this->computeTangentModuliOnQuad(args);
   }
 
   this->was_stiffness_assembled = true;
 }
 
 /* -------------------------------------------------------------------------- */
 template <Int dim>
 void MaterialElasticLinearAnisotropic<dim>::computePotentialEnergy(
     ElementType el_type) {
   AKANTU_DEBUG_ASSERT(!this->finite_deformation,
                       "finite deformation not possible in material anisotropic "
                       "(TO BE IMPLEMENTED)");
 
   auto && arguments = Material::getArguments<dim>(el_type, _not_ghost);
 
   for (auto && args :
        zip(arguments, this->potential_energy(el_type, _not_ghost))) {
     this->computePotentialEnergyOnQuad(std::get<0>(args), std::get<1>(args));
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <Int dim>
 void MaterialElasticLinearAnisotropic<dim>::computePotentialEnergyByElement(
     ElementType type, Int index, Vector<Real> & epot_on_quad_points) {
 
   auto gradu_view = make_view<dim, dim>(this->gradu(type));
   auto stress_view = make_view<dim, dim>(this->stress(type));
 
   auto nb_quadrature_points = this->fem.getNbIntegrationPoints(type);
 
   auto gradu_it = gradu_view.begin() + index * nb_quadrature_points;
   auto gradu_end = gradu_it + nb_quadrature_points;
   auto stress_it = stress_view.begin() + index * nb_quadrature_points;
   auto stress_end = stress_it + nb_quadrature_points;
 
   auto epot_quad = epot_on_quad_points.begin();
 
   Matrix<Real> grad_u(dim, dim);
 
   for (auto data : zip(tuple::get<"grad_u"_h>() = range(gradu_it, gradu_end),
                    tuple::get<"sigma"_h>() = range(stress_it, stress_end),
                    tuple::get<"Epot"_h>() = epot_on_quad_points)) {
     this->computePotentialEnergyOnQuad(data, tuple::get<"Epot"_h>(data));
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <Int dim>
 Real MaterialElasticLinearAnisotropic<dim>::getCelerity(
     const Element & /*element*/) const {
   return std::sqrt(this->eigC(0) / rho);
 }
 
 /* -------------------------------------------------------------------------- */
 
 INSTANTIATE_MATERIAL(elastic_anisotropic, MaterialElasticLinearAnisotropic);
 
 } // namespace akantu
diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc
index 4a127609b..09da84bae 100644
--- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc
+++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc
@@ -1,113 +1,113 @@
 /**
  * @file   test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc
  *
  * @author Fabian Barras <fabian.barras@epfl.ch>
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  *
  * @date creation: Sun Oct 19 2014
  * @date last modification:  Thu Feb 25 2021
  *
  * @brief  Computation of the analytical exemple 1.1 in the TGC vol 6
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 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 "sparse_matrix.hh"
 #include "test_structural_mechanics_model_fixture.hh"
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 
 using namespace akantu;
 
 /* -------------------------------------------------------------------------- */
 class TestStructDKT18 : public TestStructuralFixture<
                             element_type_t<_discrete_kirchhoff_triangle_18>> {
   using parent =
       TestStructuralFixture<element_type_t<_discrete_kirchhoff_triangle_18>>;
 
 public:
   void addMaterials() override {
     mat.E = 1;
     mat.t = 1;
     mat.nu = 0.3;
 
     this->model->addMaterial(mat);
   }
 
   void assignMaterials() override {
     auto & materials = this->model->getElementMaterial(parent::type);
     std::fill(materials.begin(), materials.end(), 0);
   }
 
   void setDirichletBCs() override {
     this->model->getBlockedDOFs().set(true);
     auto center_node = this->model->getBlockedDOFs().end(parent::ndof) - 1;
-    *center_node = Vector<bool>{false, false, false, false, false, true};
+    *center_node << false, false, false, false, false, true;
 
     this->model->getDisplacement().zero();
     auto disp = ++this->model->getDisplacement().begin(parent::ndof);
 
     // Displacement field from Batoz Vol. 2 p. 392
     // with theta to beta correction from discrete Kirchhoff condition
     // see also the master thesis of Michael Lozano
 
     // clang-format off
 
     // This displacement field tests membrane and bending modes
-    *disp = Vector<Real>{40, 20, -800 , -20, 40, 0}; ++disp;
-    *disp = Vector<Real>{50, 40, -1400, -40, 50, 0}; ++disp;
-    *disp = Vector<Real>{10, 20, -200 , -20, 10, 0}; ++disp;
+    *disp << 40, 20, -800 , -20, 40, 0; ++disp;
+    *disp << 50, 40, -1400, -40, 50, 0; ++disp;
+    *disp << 10, 20, -200 , -20, 10, 0; ++disp;
 
     // This displacement tests the bending mode
     // *disp = {0, 0, -800 , -20, 40, 0}; ++disp;
     // *disp = {0, 0, -1400, -40, 50, 0}; ++disp;
     // *disp = {0, 0, -200 , -20, 10, 0}; ++disp;
 
     // This displacement tests the membrane mode
     // *disp = {40, 20, 0, 0, 0, 0}; ++disp;
     // *disp = {50, 40, 0, 0, 0, 0}; ++disp;
     // *disp = {10, 20, 0, 0, 0, 0}; ++disp;
 
     // clang-format on
   }
 
   void setNeumannBCs() override {}
 
 protected:
   StructuralMaterial mat;
 };
 
 /* -------------------------------------------------------------------------- */
 
 // Batoz Vol 2. patch test, ISBN 2-86601-259-3
 TEST_F(TestStructDKT18, TestDisplacements) {
   model->solveStep();
   Vector<Real> solution = Vector<Real>{22.5, 22.5, -337.5, -22.5, 22.5, 0};
   auto nb_nodes = this->model->getDisplacement().size();
 
   Vector<Real> center_node_disp =
       model->getDisplacement().begin(solution.size())[nb_nodes - 1];
 
   auto error = solution - center_node_disp;
 
   EXPECT_NEAR(error.norm(), 0., 1e-12);
 }