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 98b0cb329..ed3baeb2f 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,313 +1,879 @@
 /* -------------------------------------------------------------------------- */
 #include "material_elastic.hh"
 #include "material_elastic_orthotropic.hh"
 #include "solid_mechanics_model.hh"
 #include "test_material_fixtures.hh"
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 using types =
     ::testing::Types<Traits<MaterialElastic, 1>, Traits<MaterialElastic, 2>,
                      Traits<MaterialElastic, 3>,
 
                      Traits<MaterialElasticOrthotropic, 1>,
                      Traits<MaterialElasticOrthotropic, 2>,
                      Traits<MaterialElasticOrthotropic, 3>,
 
                      Traits<MaterialElasticLinearAnisotropic, 1>,
                      Traits<MaterialElasticLinearAnisotropic, 2>,
                      Traits<MaterialElasticLinearAnisotropic, 3>>;
 
 /* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<3>>::testPushWaveSpeed() {
+template <> void FriendMaterial<MaterialElastic<1>>::testComputeStress() {
+  Real E = 3.;
+  setParam("E", E);
+
+  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;
+  ASSERT_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;
+  ASSERT_NEAR(epot, solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <>
+void FriendMaterial<MaterialElastic<1>>::testComputeTangentModuli() {
+  Real E = 2;
+  setParam("E", E);
+  Matrix<Real> tangent(1, 1);
+  this->computeTangentModuliOnQuad(tangent);
+  ASSERT_NEAR(tangent(0, 0), E, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<1>>::testPushWaveSpeed() {
+  Real E = 3., rho = 2.;
+  setParam("E", E);
+  setParam("rho", rho);
+
+  auto wave_speed = this->getPushWaveSpeed(Element());
+  auto solution = std::sqrt(E / rho);
+  ASSERT_NEAR(wave_speed, solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<1>>::testShearWaveSpeed() {
+  ASSERT_THROW(this->getShearWaveSpeed(Element()), debug::Exception);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<2>>::testComputeStress() {
+  Real E = 1.;
+  Real nu = .3;
+  Real sigma_th = 0.3; // thermal stress
+  setParam("E", E);
+  setParam("nu", nu);
+
+  Matrix<Real> rotation_matrix = getRandomRotation2d();
+
+  auto grad_u = this->getDeviatoricStrain(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);
+
+  Matrix<Real> identity(2, 2);
+  identity.eye();
+
+  Matrix<Real> sigma_expected =
+    0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity;
+
+  auto diff = sigma - sigma_expected;
+  Real stress_error = diff.norm<L_inf>();
+  ASSERT_NEAR(stress_error, 0., 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+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);
+  ASSERT_NEAR(epot, solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <>
+void FriendMaterial<MaterialElastic<2>>::testComputeTangentModuli() {
+  Real E = 1.;
+  Real nu = .3;
+  setParam("E", E);
+  setParam("nu", nu);
+  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>();
+  ASSERT_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>();
+  ASSERT_NEAR(tangent_error, 0, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<2>>::testPushWaveSpeed() {
   Real E = 1.;
   Real nu = .3;
   Real rho = 2;
   setParam("E", E);
   setParam("nu", nu);
   setParam("rho", rho);
   auto wave_speed = this->getPushWaveSpeed(Element());
 
   Real K = E / (3 * (1 - 2 * nu));
   Real mu = E / (2 * (1 + nu));
   Real sol = std::sqrt((K + 4. / 3 * mu) / rho);
 
   ASSERT_NEAR(wave_speed, sol, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<3>>::testShearWaveSpeed() {
+template <> void FriendMaterial<MaterialElastic<2>>::testShearWaveSpeed() {
   Real E = 1.;
   Real nu = .3;
   Real rho = 2;
   setParam("E", E);
   setParam("nu", nu);
   setParam("rho", rho);
   auto wave_speed = this->getShearWaveSpeed(Element());
 
   Real mu = E / (2 * (1 + nu));
   Real sol = std::sqrt(mu / rho);
 
   ASSERT_NEAR(wave_speed, sol, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> void FriendMaterial<MaterialElastic<3>>::testComputeStress() {
   Real E = 1.;
   Real nu = .3;
   Real sigma_th = 0.3; // thermal stress
   setParam("E", E);
   setParam("nu", nu);
 
   Matrix<Real> rotation_matrix = getRandomRotation3d();
 
   auto grad_u = this->getDeviatoricStrain(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> sigma_expected =
       0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity;
 
   auto diff = sigma - sigma_expected;
   Real stress_error = diff.norm<L_inf>();
   ASSERT_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);
+  ASSERT_NEAR(epot, solution, 1e-14);
+}
+
 /* -------------------------------------------------------------------------- */
 template <>
 void FriendMaterial<MaterialElastic<3>>::testComputeTangentModuli() {
   Real E = 1.;
   Real nu = .3;
   setParam("E", E);
   setParam("nu", nu);
   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>();
   ASSERT_NEAR(tangent_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);
-  ASSERT_NEAR(epot, solution, 1e-14);
-}
-
-/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<1>>::testComputeStress() {
-  Real E = 3.;
-  setParam("E", E);
-
-  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;
-  ASSERT_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;
-  ASSERT_NEAR(epot, solution, 1e-14);
-}
-
-/* -------------------------------------------------------------------------- */
-template <>
-void FriendMaterial<MaterialElastic<1>>::testComputeTangentModuli() {
-  Real E = 2;
-  setParam("E", E);
-  Matrix<Real> tangent(1, 1);
-  this->computeTangentModuliOnQuad(tangent);
-  ASSERT_NEAR(tangent(0, 0), E, 1e-14);
-}
-
-/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<1>>::testPushWaveSpeed() {
-  Real E = 3., rho = 2.;
-  setParam("E", E);
-  setParam("rho", rho);
-
-  auto wave_speed = this->getPushWaveSpeed(Element());
-  auto solution = std::sqrt(E / rho);
-  ASSERT_NEAR(wave_speed, solution, 1e-14);
-}
-
-/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<1>>::testShearWaveSpeed() {
-  ASSERT_THROW(this->getShearWaveSpeed(Element()), debug::Exception);
-}
-
-/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<2>>::testPushWaveSpeed() {
+template <> void FriendMaterial<MaterialElastic<3>>::testPushWaveSpeed() {
   Real E = 1.;
   Real nu = .3;
   Real rho = 2;
   setParam("E", E);
   setParam("nu", nu);
   setParam("rho", rho);
   auto wave_speed = this->getPushWaveSpeed(Element());
 
   Real K = E / (3 * (1 - 2 * nu));
   Real mu = E / (2 * (1 + nu));
   Real sol = std::sqrt((K + 4. / 3 * mu) / rho);
 
   ASSERT_NEAR(wave_speed, sol, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<2>>::testShearWaveSpeed() {
+template <> void FriendMaterial<MaterialElastic<3>>::testShearWaveSpeed() {
   Real E = 1.;
   Real nu = .3;
   Real rho = 2;
   setParam("E", E);
   setParam("nu", nu);
   setParam("rho", rho);
   auto wave_speed = this->getShearWaveSpeed(Element());
 
   Real mu = E / (2 * (1 + nu));
   Real sol = std::sqrt(mu / rho);
 
   ASSERT_NEAR(wave_speed, sol, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<2>>::testComputeStress() {
-  Real E = 1.;
-  Real nu = .3;
-  Real sigma_th = 0.3; // thermal stress
-  setParam("E", E);
-  setParam("nu", nu);
-
+template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testComputeStress() {
+  UInt Dim = 2;
+  Real E1 = 1.;
+  Real E2 = 2.;
+  Real nu12 = 0.1;
+  Real G12 = 2.;
+
+  setParamNoUpdate("E1", E1);
+  setParamNoUpdate("E2", E2);
+  setParamNoUpdate("nu12", nu12);
+  setParamNoUpdate("G12", G12);
+
+  // material frame of reference is rotate by rotation_matrix starting from
+  // canonical basis
   Matrix<Real> rotation_matrix = getRandomRotation2d();
 
-  auto grad_u = this->getDeviatoricStrain(1.).block(0, 0, 2, 2);
+  // 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)
+  Vector<Real> v1(Dim);
+  Vector<Real> v2(Dim);
+  for (UInt i = 0; i < Dim; ++i) {
+    v1[i] = rotation_matrix(i, 0);
+    v2[i] = rotation_matrix(i, 1);
+  }
+
+  setParamNoUpdate("n1", v1.normalize());
+  setParamNoUpdate("n2", v2.normalize());
+
+  // set internal Cijkl matrix expressed in the canonical frame of reference
+  this->updateInternalParameters();
 
-  auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix);
+  // gradient in material frame of reference
+  auto grad_u = ( this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(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, sigma_th);
+  this->computeStressOnQuad(grad_u_rot, sigma_rot);
 
-  auto sigma = this->reverseRotation(sigma_rot, rotation_matrix);
+// stress in the material reference (we need to apply the rotation)
+  auto sigma = this->applyRotation(sigma_rot, rotation_matrix);
 
-  Matrix<Real> identity(2, 2);
-  identity.eye();
+// 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> sigma_expected =
-      0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity;
+  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>();
-  ASSERT_NEAR(stress_error, 0., 1e-14);
+  ASSERT_NEAR(stress_error, 0., 1e-13);
 }
 
 /* -------------------------------------------------------------------------- */
-template <>
-void FriendMaterial<MaterialElastic<2>>::testComputeTangentModuli() {
-  Real E = 1.;
-  Real nu = .3;
-  setParam("E", E);
-  setParam("nu", nu);
-  Matrix<Real> tangent(3, 3);
+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);
+  ASSERT_NEAR(epot, solution, 1e-14);
+}
 
-  /* 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));
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testComputeTangentModuli() {
+
+  // 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.;
+
+  setParamNoUpdate("n1", n1);
+  setParamNoUpdate("n2", n2);
+  setParamNoUpdate("E1", E1);
+  setParamNoUpdate("E2", E2);
+  setParamNoUpdate("nu12", nu12);
+  setParamNoUpdate("G12", G12);
+
+  // set internal Cijkl matrix expressed in the canonical frame of reference
+  this->updateInternalParameters();
 
+// 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 - solution).norm<L_2>();
+
+  Real tangent_error = (tangent - C_expected).norm<L_2>();
   ASSERT_NEAR(tangent_error, 0, 1e-14);
+}
 
-  /* Plane Stress */
-  this->plane_stress = true;
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testComputeStress() {
+  UInt Dim = 3;
+  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.;
+
+  setParamNoUpdate("E1", E1);
+  setParamNoUpdate("E2", E2);
+  setParamNoUpdate("E3", E3);
+  setParamNoUpdate("nu12", nu12);
+  setParamNoUpdate("nu13", nu13);
+  setParamNoUpdate("nu23", nu23);
+  setParamNoUpdate("G12", G12);
+  setParamNoUpdate("G13", G13);
+  setParamNoUpdate("G23", G23);
+
+  // material frame of reference is rotate by rotation_matrix starting from
+  // canonical basis
+  Matrix<Real> rotation_matrix = getRandomRotation3d();
+
+  // 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)
+  Vector<Real> v1(Dim);
+  Vector<Real> v2(Dim);
+  Vector<Real> v3(Dim);
+  for (UInt i = 0; i < Dim; ++i) {
+    v1[i] = rotation_matrix(i, 0);
+    v2[i] = rotation_matrix(i, 1);
+    v3[i] = rotation_matrix(i, 2);
+  }
+
+  setParamNoUpdate("n1", v1.normalize());
+  setParamNoUpdate("n2", v2.normalize());
+  setParamNoUpdate("n3", v3.normalize());
+
+  // set internal Cijkl matrix expressed in the canonical frame of reference
   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);
 
+  // gradient in material frame of reference
+  auto grad_u = this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.);
+
+  // 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>();
+  ASSERT_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);
+  ASSERT_NEAR(epot, solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testComputeTangentModuli() {
+
+  // Note: for this test material and canonical basis coincide
+  UInt Dim = 3;
+  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.;
+
+  setParamNoUpdate("n1", n1);
+  setParamNoUpdate("n2", n2);
+  setParamNoUpdate("n3", n3);
+  setParamNoUpdate("E1", E1);
+  setParamNoUpdate("E2", E2);
+  setParamNoUpdate("E3", E3);
+  setParamNoUpdate("nu12", nu12);
+  setParamNoUpdate("nu13", nu13);
+  setParamNoUpdate("nu23", nu23);
+  setParamNoUpdate("G12", G12);
+  setParamNoUpdate("G13", G13);
+  setParamNoUpdate("G23", G23);
+
+  // set internal Cijkl matrix expressed in the canonical frame of reference
+  this->updateInternalParameters();
+
+// 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);
-  tangent_error = (tangent - solution).norm<L_2>();
+
+  Real tangent_error = (tangent - C_expected).norm<L_2>();
   ASSERT_NEAR(tangent_error, 0, 1e-14);
 }
 
 /* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<2>>::testEnergyDensity() {
+template <> void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testComputeStress() {
+  UInt Dim = 2;
+  Matrix<Real> C = {
+    {1.0, 0.3, 0.4},
+    {0.3, 2.0, 0.1},
+    {0.4, 0.1, 1.5},
+  };
+
+  setParamNoUpdate("C11", C(0,0));
+  setParamNoUpdate("C12", C(0,1));
+  setParamNoUpdate("C13", C(0,2));
+  setParamNoUpdate("C22", C(1,1));
+  setParamNoUpdate("C23", C(1,2));
+  setParamNoUpdate("C33", C(2,2));
+
+  // material frame of reference is rotate by rotation_matrix starting from
+  // canonical basis
+  Matrix<Real> rotation_matrix = getRandomRotation2d();
+
+  // 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)
+  Vector<Real> v1(Dim);
+  Vector<Real> v2(Dim);
+  for (UInt i = 0; i < Dim; ++i) {
+    v1[i] = rotation_matrix(i, 0);
+    v2[i] = rotation_matrix(i, 1);
+  }
+
+  setParamNoUpdate("n1", v1.normalize());
+  setParamNoUpdate("n2", v2.normalize());
+
+  // set internal Cijkl matrix expressed in the canonical frame of reference
+  this->updateInternalParameters();
+
+  // gradient in material frame of reference
+  auto grad_u = (this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(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(Dim, Dim);
+  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>();
+  ASSERT_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);
   ASSERT_NEAR(epot, solution, 1e-14);
 }
 
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testComputeTangentModuli() {
+
+  // Note: for this test material and canonical basis coincide
+  Matrix<Real> C = {
+    {1.0, 0.3, 0.4},
+    {0.3, 2.0, 0.1},
+    {0.4, 0.1, 1.5},
+  };
+
+  setParamNoUpdate("C11", C(0,0));
+  setParamNoUpdate("C12", C(0,1));
+  setParamNoUpdate("C13", C(0,2));
+  setParamNoUpdate("C22", C(1,1));
+  setParamNoUpdate("C23", C(1,2));
+  setParamNoUpdate("C33", C(2,2));
+
+  // set internal Cijkl matrix expressed in the canonical frame of reference
+  this->updateInternalParameters();
+
+  Matrix<Real> tangent(3,3);
+  this->computeTangentModuliOnQuad(tangent);
+
+  Real tangent_error = (tangent - C).norm<L_2>();
+  ASSERT_NEAR(tangent_error, 0, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testComputeStress() {
+  UInt Dim = 3;
+  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},
+  };
+
+  setParamNoUpdate("C11", C(0,0));
+  setParamNoUpdate("C12", C(0,1));
+  setParamNoUpdate("C13", C(0,2));
+  setParamNoUpdate("C14", C(0,3));
+  setParamNoUpdate("C15", C(0,4));
+  setParamNoUpdate("C16", C(0,5));
+  setParamNoUpdate("C22", C(1,1));
+  setParamNoUpdate("C23", C(1,2));
+  setParamNoUpdate("C24", C(1,3));
+  setParamNoUpdate("C25", C(1,4));
+  setParamNoUpdate("C26", C(1,5));
+  setParamNoUpdate("C33", C(2,2));
+  setParamNoUpdate("C34", C(2,3));
+  setParamNoUpdate("C35", C(2,4));
+  setParamNoUpdate("C36", C(2,5));
+  setParamNoUpdate("C44", C(3,3));
+  setParamNoUpdate("C45", C(3,4));
+  setParamNoUpdate("C46", C(3,5));
+  setParamNoUpdate("C55", C(4,4));
+  setParamNoUpdate("C56", C(4,5));
+  setParamNoUpdate("C66", C(5,5));
+
+  // material frame of reference is rotate by rotation_matrix starting from
+  // canonical basis
+  Matrix<Real> rotation_matrix = getRandomRotation3d();
+
+  // 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)
+  Vector<Real> v1(Dim);
+  Vector<Real> v2(Dim);
+  Vector<Real> v3(Dim);
+  for (UInt i = 0; i < Dim; ++i) {
+    v1[i] = rotation_matrix(i, 0);
+    v2[i] = rotation_matrix(i, 1);
+    v3[i] = rotation_matrix(i, 2);
+  }
+
+  setParamNoUpdate("n1", v1.normalize());
+  setParamNoUpdate("n2", v2.normalize());
+  setParamNoUpdate("n3", v3.normalize());
+
+  // set internal Cijkl matrix expressed in the canonical frame of reference
+  this->updateInternalParameters();
+
+  // gradient in material frame of reference
+  auto grad_u = this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.);
+
+  // 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(Dim, Dim);
+  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>();
+  ASSERT_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);
+  ASSERT_NEAR(epot, solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testComputeTangentModuli() {
+
+  // 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},
+  };
+
+  setParamNoUpdate("C11", C(0,0));
+  setParamNoUpdate("C12", C(0,1));
+  setParamNoUpdate("C13", C(0,2));
+  setParamNoUpdate("C14", C(0,3));
+  setParamNoUpdate("C15", C(0,4));
+  setParamNoUpdate("C16", C(0,5));
+  setParamNoUpdate("C22", C(1,1));
+  setParamNoUpdate("C23", C(1,2));
+  setParamNoUpdate("C24", C(1,3));
+  setParamNoUpdate("C25", C(1,4));
+  setParamNoUpdate("C26", C(1,5));
+  setParamNoUpdate("C33", C(2,2));
+  setParamNoUpdate("C34", C(2,3));
+  setParamNoUpdate("C35", C(2,4));
+  setParamNoUpdate("C36", C(2,5));
+  setParamNoUpdate("C44", C(3,3));
+  setParamNoUpdate("C45", C(3,4));
+  setParamNoUpdate("C46", C(3,5));
+  setParamNoUpdate("C55", C(4,4));
+  setParamNoUpdate("C56", C(4,5));
+  setParamNoUpdate("C66", C(5,5));
+
+  // set internal Cijkl matrix expressed in the canonical frame of reference
+  this->updateInternalParameters();
+
+  Matrix<Real> tangent(6,6);
+  this->computeTangentModuliOnQuad(tangent);
+
+  Real tangent_error = (tangent - C).norm<L_2>();
+  ASSERT_NEAR(tangent_error, 0, 1e-14);
+}
+
 /* -------------------------------------------------------------------------- */
 namespace {
 
 template <typename T>
 class TestElasticMaterialFixture : public ::TestMaterialFixture<T> {};
 
 TYPED_TEST_CASE(TestElasticMaterialFixture, types);
 
 TYPED_TEST(TestElasticMaterialFixture, ElasticComputeStress) {
   this->material->testComputeStress();
 }
 
 TYPED_TEST(TestElasticMaterialFixture, ElasticEnergyDensity) {
   this->material->testEnergyDensity();
 }
 
 TYPED_TEST(TestElasticMaterialFixture, ElasticComputeTangentModuli) {
   this->material->testComputeTangentModuli();
 }
 
 TYPED_TEST(TestElasticMaterialFixture, ElasticComputePushWaveSpeed) {
   this->material->testPushWaveSpeed();
 }
 
 TYPED_TEST(TestElasticMaterialFixture, ElasticComputeShearWaveSpeed) {
   this->material->testShearWaveSpeed();
 }
 } // namespace
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh b/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh
index 4d4129a36..c0a846674 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh
@@ -1,180 +1,181 @@
 /* -------------------------------------------------------------------------- */
 #include "material_elastic.hh"
 #include "solid_mechanics_model.hh"
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 #include <random>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 #define TO_IMPLEMENT AKANTU_EXCEPTION("TEST TO IMPLEMENT")
 
 using namespace akantu;
 
 /* -------------------------------------------------------------------------- */
 template <typename T> class FriendMaterial : public T {
 public:
   ~FriendMaterial() = default;
   
   virtual void testComputeStress() { TO_IMPLEMENT; };
   virtual void testComputeTangentModuli() { TO_IMPLEMENT; };
   virtual void testEnergyDensity() { TO_IMPLEMENT; };
   virtual void testPushWaveSpeed() { TO_IMPLEMENT; }
   virtual void testShearWaveSpeed() { TO_IMPLEMENT; }
 
   static inline Matrix<Real> getDeviatoricStrain(Real intensity);
 
   static inline Matrix<Real> getHydrostaticStrain(Real intensity);
 
   static inline const Matrix<Real>
   reverseRotation(Matrix<Real> mat, Matrix<Real> rotation_matrix) {
     Matrix<Real> R = rotation_matrix;
     Matrix<Real> m2 = mat * R;
     Matrix<Real> m1 = R.transpose();
     return m1 * m2;
   };
 
   static inline const Matrix<Real>
   applyRotation(Matrix<Real> mat, const Matrix<Real> rotation_matrix) {
     Matrix<Real> R = rotation_matrix;
     Matrix<Real> m2 = mat * R.transpose();
     Matrix<Real> m1 = R;
     return m1 * m2;
   };
 
   static inline Matrix<Real> getRandomRotation3d();
   static inline Matrix<Real> getRandomRotation2d();
 
   using T::T;
 };
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 Matrix<Real> FriendMaterial<T>::getDeviatoricStrain(Real intensity) {
   Matrix<Real> dev = {{0, 1, 2}, {1, 0, 3}, {2, 3, 0}};
   dev *= intensity;
   return dev;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 Matrix<Real> FriendMaterial<T>::getHydrostaticStrain(Real intensity) {
   Matrix<Real> dev = {{1, 0, 0}, {0, 2, 0}, {0, 0, 3}};
   dev *= intensity;
   return dev;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> Matrix<Real> FriendMaterial<T>::getRandomRotation3d() {
   constexpr UInt Dim = 3;
   // Seed with a real random value, if available
   std::random_device rd;
   std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd()
   std::uniform_real_distribution<> dis;
 
   Vector<Real> v1(Dim);
   Vector<Real> v2(Dim);
   Vector<Real> v3(Dim);
 
   for (UInt i = 0; i < Dim; ++i) {
     v1(i) = dis(gen);
     v2(i) = dis(gen);
   }
 
   v1.normalize();
   v2.normalize();
 
   auto d = v1.dot(v2);
   v2 -= v1 * d;
   v2.normalize();
 
   d = v1.dot(v2);
   v2 -= v1 * d;
   v2.normalize();
 
   v3.crossProduct(v1, v2);
 
   Vector<Real> test_axis(3);
   test_axis.crossProduct(v1, v2);
   test_axis -= v3;
   if (test_axis.norm() > 8 * std::numeric_limits<Real>::epsilon()) {
     AKANTU_DEBUG_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() << ".");
   }
 
   Matrix<Real> mat(Dim, Dim);
   for (UInt i = 0; i < Dim; ++i) {
     mat(0, i) = v1[i];
     mat(1, i) = v2[i];
     mat(2, i) = v3[i];
   }
 
   return mat;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> Matrix<Real> FriendMaterial<T>::getRandomRotation2d() {
   constexpr UInt Dim = 2;
   // Seed with a real random value, if available
   std::random_device rd;
   std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd()
   std::uniform_real_distribution<> dis;
 
   // v1 and v2 are such that they form a right-hand basis with prescribed v3,
-  // it's need (at least) for the 2d lin. elastic orthotropic material
+  // it's need (at least) for 2d linear elastic anisotropic and (orthotropic)
+  // materials to rotate the tangent stiffness
   
   Vector<Real> v1(3);
   Vector<Real> v2(3);
   Vector<Real> v3 = {0,0,1};
 
   for (UInt i = 0; i < Dim; ++i) {
     v1(i) = dis(gen);
     // v2(i) = dis(gen);
   }
 
   v1.normalize();
   v2.crossProduct(v3,v1);
 
   Matrix<Real> mat(Dim, Dim);
   for (UInt i = 0; i < Dim; ++i) {
     mat(0, i) = v1[i];
     mat(1, i) = v2[i];
   }
 
   return mat;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> class TestMaterialFixture : public ::testing::Test {
 public:
   using mat_class = typename T::mat_class;
 
   void SetUp() override {
     constexpr auto spatial_dimension = T::Dim;
     mesh = std::make_unique<Mesh>(spatial_dimension);
     model = std::make_unique<SolidMechanicsModel>(*mesh);
     material = std::make_unique<friend_class>(*model);
   }
 
   void TearDown() override {
     material.reset(nullptr);
     model.reset(nullptr);
     mesh.reset(nullptr);
   }
 
   using friend_class = FriendMaterial<mat_class>;
 
 protected:
   std::unique_ptr<Mesh> mesh;
   std::unique_ptr<SolidMechanicsModel> model;
   std::unique_ptr<friend_class> material;
 };
 
 /* -------------------------------------------------------------------------- */
 template <template <UInt> class T, UInt _Dim> struct Traits {
   using mat_class = T<_Dim>;
   static constexpr UInt Dim = _Dim;
 };
 
 /* -------------------------------------------------------------------------- */