diff --git a/elements/src/TimoshenkoBeamElement.cc b/elements/src/TimoshenkoBeamElement.cc
index bdfdcde5746ba8d4a4058a8aa3d39e10d0d281e4..b37e16a945f48fc1f990cc7afd46cad08ccf4da6 100644
--- a/elements/src/TimoshenkoBeamElement.cc
+++ b/elements/src/TimoshenkoBeamElement.cc
@@ -12,4 +12,204 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
-#include "ae108/elements/TimoshenkoBeamElement.h"
\ No newline at end of file
+#include "ae108/elements/TimoshenkoBeamElement.h"
+
+namespace ae108 {
+namespace elements {
+namespace {
+
+/**
+ * @brief Computes the stiffness matrix of a reference beam with the given
+ * properties.
+ */
+template <std::size_t Dimension_>
+Eigen::Matrix<double, Dimension_ *(Dimension_ + 1),
+              Dimension_ *(Dimension_ + 1), Eigen::RowMajor>
+stiffness_matrix(const Properties<double, Dimension_> &properties,
+                 const double length) noexcept;
+
+// refer to Cook et. al (2002), "Concepts and applications of Finite Element
+// Analysis", 4th ed., p.27
+template <>
+Eigen::Matrix<double, 12, 12, Eigen::RowMajor>
+stiffness_matrix<3>(const Properties<double, 3> &properties,
+                    const double length) noexcept {
+  const auto L = length;
+  const auto A = properties.area;
+  const auto E = properties.young_modulus;
+  const auto G = properties.shear_modulus;
+  const auto I_z = properties.area_moment_z;
+  const auto k_y = properties.shear_correction_factor_y;
+  const auto I_y = properties.area_moment_y;
+  const auto J_x = properties.polar_moment_x;
+  const auto k_z = properties.shear_correction_factor_z;
+
+  const auto phi = [&](const double I, const double k) {
+    return 12. * E * I * k / A / G / L / L;
+  };
+  const double phi_y = phi(I_z, k_y);
+  const double phi_z = phi(I_y, k_z);
+
+  const auto f1 = [&](const double I, const double phi) {
+    return 12. * E * I / (1. + phi) / L / L / L;
+  };
+  const auto Y1 = f1(I_z, phi_y);
+  const auto Z1 = f1(I_y, phi_z);
+
+  const auto f2 = [&](const double I, const double phi) {
+    return 6. * E * I / (1. + phi) / L / L;
+  };
+  const auto Y2 = f2(I_z, phi_y);
+  const auto Z2 = f2(I_y, phi_z);
+
+  const auto f3 = [&](const double I, const double phi) {
+    return (4. + phi) * E * I / (1. + phi) / L;
+  };
+  const auto Y3 = f3(I_z, phi_y);
+  const auto Z3 = f3(I_y, phi_z);
+
+  const auto f4 = [&](const double I, const double phi) {
+    return (2. - phi) * E * I / (1. + phi) / L;
+  };
+  const auto Y4 = f4(I_z, phi_y);
+  const auto Z4 = f4(I_y, phi_z);
+
+  const auto X = A * E / L;
+  const auto S = G * J_x / L;
+  const auto _ = 0.;
+
+  // clang-format off
+  const tensor::Tensor<double, 12, 12> matrix = {{
+      {{  X,   _,   _,   _,   _,   _,  -X,   _,   _,   _,   _,   _}},
+      {{  _,  Y1,   _,   _,   _,  Y2,   _, -Y1,   _,   _,   _,  Y2}},
+      {{  _,   _,  Z1,   _, -Z2,   _,   _,   _, -Z1,   _, -Z2,   _}},
+      {{  _,   _,   _,   S,   _,   _,   _,   _,   _,  -S,   _,   _}},
+      {{  _,   _, -Z2,   _,  Z3,   _,   _,   _,  Z2,   _,  Z4,   _}},
+      {{  _,  Y2,   _,   _,   _,  Y3,   _, -Y2,   _,   _,   _,  Y4}},
+      {{ -X,   _,   _,   _,   _,   _,   X,   _,   _,   _,   _,   _}},
+      {{  _, -Y1,   _,   _,   _, -Y2,   _,  Y1,   _,   _,   _, -Y2}},
+      {{  _,   _, -Z1,   _,  Z2,   _,   _,   _,  Z1,   _,  Z2,   _}},
+      {{  _,   _,   _,  -S,   _,   _,   _,   _,   _,   S,   _,   _}},
+      {{  _,   _, -Z2,   _,  Z4,   _,   _,   _,  Z2,   _,  Z3,   _}},
+      {{  _,  Y2,   _,   _,   _,  Y4,   _, -Y2,   _,   _,   _,  Y3}},
+  }};
+  // clang-format on
+
+  return tensor::as_matrix_of_rows(&matrix);
+}
+
+// refer to Cook et. al (2002), "Concepts and applications of Finite Element
+// Analysis", 4th ed., p.26
+template <>
+Eigen::Matrix<double, 6, 6, Eigen::RowMajor>
+stiffness_matrix<2>(const Properties<double, 2> &properties,
+                    const double length) noexcept {
+  const auto L = length;
+  const auto A = properties.area;
+  const auto E = properties.young_modulus;
+  const auto G = properties.shear_modulus;
+  const auto I_z = properties.area_moment_z;
+  const auto k_y = properties.shear_correction_factor_y;
+
+  const double phi_y = 12 * E * I_z * k_y / A / G / L / L;
+
+  const auto X = A * E / L;
+  const auto Y1 = 12 * E * I_z / (1 + phi_y) / L / L / L;
+  const auto Y2 = 6 * E * I_z / (1 + phi_y) / L / L;
+  const auto Y3 = (4 + phi_y) * E * I_z / (1 + phi_y) / L;
+  const auto Y4 = (2 - phi_y) * E * I_z / (1 + phi_y) / L;
+
+  const auto _ = 0.;
+
+  // clang-format off
+  const tensor::Tensor<double, 6, 6> matrix = {{
+      {{  X,   _,   _,  -X,   _,   _}},
+      {{  _,  Y1,  Y2,   _, -Y1,  Y2}},
+      {{  _,  Y2,  Y3,   _, -Y2,  Y4}},
+      {{ -X,   _,   _,   X,   _,   _}},
+      {{  _, -Y1, -Y2,   _,  Y1, -Y2}},
+      {{  _,  Y2,  Y4,   _, -Y2,  Y3}},
+  }};
+  // clang-format on
+
+  return tensor::as_matrix_of_rows(&matrix);
+}
+
+template <std::size_t Dimension_>
+Eigen::Matrix<double, Dimension_ *(Dimension_ + 1),
+              Dimension_ *(Dimension_ + 1), Eigen::RowMajor>
+rotation_matrix(const tensor::Tensor<double, Dimension_> &orientation) noexcept;
+
+// refer to Cook et. al (2002), "Concepts and applications of Finite Element
+// Analysis", 4th ed., p.32
+template <>
+Eigen::Matrix<double, 12, 12, Eigen::RowMajor>
+rotation_matrix<3>(const tensor::Tensor<double, 3> &orientation) noexcept {
+  // rotation that maps the normalized orientation vector to (1, 0, 0)
+  const auto Lambda = Eigen::Quaternion<double>()
+                          .FromTwoVectors(tensor::as_vector(&orientation),
+                                          Eigen::Vector3d::UnitX())
+                          .normalized()
+                          .toRotationMatrix();
+
+  auto result = Eigen::Matrix<double, 12, 12, Eigen::RowMajor>::Zero().eval();
+
+  result.block(0, 0, 3, 3) = Lambda;
+  result.block(3, 3, 3, 3) = Lambda;
+  result.block(6, 6, 3, 3) = Lambda;
+  result.block(9, 9, 3, 3) = Lambda;
+
+  return result;
+}
+
+// refer to Cook et. al (2002), "Concepts and applications of Finite Element
+// Analysis", 4th ed., p.31
+template <>
+Eigen::Matrix<double, 6, 6, Eigen::RowMajor>
+rotation_matrix<2>(const tensor::Tensor<double, 2> &orientation) noexcept {
+  const auto normalized =
+      tensor::as_vector(&orientation) / tensor::as_vector(&orientation).norm();
+
+  // rotation that maps the normalized orientation vector to (1, 0)
+  const tensor::Tensor<double, 2, 2> Lambda = {{
+      {{normalized[0], normalized[1]}},
+      {{-normalized[1], normalized[0]}},
+  }};
+
+  auto result = Eigen::Matrix<double, 6, 6, Eigen::RowMajor>::Zero().eval();
+
+  result.block(0, 0, 2, 2) = tensor::as_matrix_of_rows(&Lambda);
+  result(2, 2) = 1.;
+  result.block(3, 3, 2, 2) = tensor::as_matrix_of_rows(&Lambda);
+  result(5, 5) = 1.;
+
+  return result;
+}
+} // namespace
+
+template <std::size_t Dimension_>
+Eigen::Matrix<double, Dimension_ *(Dimension_ + 1),
+              Dimension_ *(Dimension_ + 1), Eigen::RowMajor>
+timoshenko_beam_stiffness_matrix(
+    const tensor::Tensor<double, Dimension_> &axis,
+    const Properties<double, Dimension_> &properties) noexcept {
+  const auto reference =
+      stiffness_matrix<Dimension_>(properties, tensor::as_vector(&axis).norm());
+
+  const auto rotation = rotation_matrix<Dimension_>(axis);
+
+  return rotation.transpose() * reference * rotation;
+}
+
+template Eigen::Matrix<double, 6, 6, Eigen::RowMajor>
+timoshenko_beam_stiffness_matrix(
+    const tensor::Tensor<double, 2> &axis,
+    const Properties<double, 2> &properties) noexcept;
+
+template Eigen::Matrix<double, 12, 12, Eigen::RowMajor>
+timoshenko_beam_stiffness_matrix(
+    const tensor::Tensor<double, 3> &axis,
+    const Properties<double, 3> &properties) noexcept;
+
+} // namespace elements
+} // namespace ae108
\ No newline at end of file
diff --git a/elements/src/include/ae108/elements/TimoshenkoBeamElement.h b/elements/src/include/ae108/elements/TimoshenkoBeamElement.h
index 3a1b7e1f0de270c5748f3257e522b795410454e0..5d89d84c0ebf7a0805cbd01f371e2113f3ab0fcf 100644
--- a/elements/src/include/ae108/elements/TimoshenkoBeamElement.h
+++ b/elements/src/include/ae108/elements/TimoshenkoBeamElement.h
@@ -24,7 +24,6 @@
 
 namespace ae108 {
 namespace elements {
-namespace timoshenko {
 
 template <class ValueType_, std::size_t Dimension_> struct Properties;
 
@@ -58,154 +57,17 @@ template <class ValueType_> struct Properties<ValueType_, 2> {
 };
 
 /**
- * @brief Computes the stiffness matrix of a reference beam with the given
- * properties.
+ * @brief Computes the stiffness matrix for a Timoshenko beam with the given
+ * axis and the given properties.
+ * @tparam Dimension_ The dimenion of the physical space. Only dimensions 2 and
+ * 3 are supported.
  */
-template <class ValueType_, std::size_t Dimension_>
-Eigen::Matrix<ValueType_, Dimension_ *(Dimension_ + 1),
-              Dimension_ *(Dimension_ + 1), Eigen::RowMajor>
-stiffness_matrix(const Properties<ValueType_, Dimension_> &properties,
-                 const ValueType_ length);
-
-// refer to Cook et. al (2002), "Concepts and applications of Finite Element
-// Analysis", 4th ed., p.27
-template <>
-inline Eigen::Matrix<double, 12, 12, Eigen::RowMajor>
-stiffness_matrix<double, 3>(const Properties<double, 3> &properties,
-                            const double length) {
-  const auto L = length;
-  const auto A = properties.area;
-  const auto E = properties.young_modulus;
-  const auto G = properties.shear_modulus;
-  const auto I_z = properties.area_moment_z;
-  const auto k_y = properties.shear_correction_factor_y;
-  const auto I_y = properties.area_moment_y;
-  const auto J_x = properties.polar_moment_x;
-  const auto k_z = properties.shear_correction_factor_z;
-
-  const double phi_y = 12 * E * I_z * k_y / A / G / L / L;
-  const double phi_z = 12 * E * I_y * k_z / A / G / L / L;
-
-  const auto X = A * E / L;
-  const auto Y1 = 12 * E * I_z / (1 + phi_y) / L / L / L;
-  const auto Y2 = 6 * E * I_z / (1 + phi_y) / L / L;
-  const auto Y3 = (4 + phi_y) * E * I_z / (1 + phi_y) / L;
-  const auto Y4 = (2 - phi_y) * E * I_z / (1 + phi_y) / L;
-  const auto Z1 = 12 * E * I_y / (1 + phi_z) / L / L / L;
-  const auto Z2 = 6 * E * I_y / (1 + phi_z) / L / L;
-  const auto Z3 = (4 + phi_z) * E * I_y / (1 + phi_z) / L;
-  const auto Z4 = (2 - phi_z) * E * I_y / (1 + phi_z) / L;
-  const auto S = G * J_x / L;
-
-  const auto _ = 0.;
-
-  // clang-format off
-  const tensor::Tensor<double, 12, 12> matrix = {{
-      {{  X,   _,   _,   _,   _,   _,  -X,   _,   _,   _,   _,   _}},
-      {{  _,  Y1,   _,   _,   _,  Y2,   _, -Y1,   _,   _,   _,  Y2}},
-      {{  _,   _,  Z1,   _, -Z2,   _,   _,   _, -Z1,   _, -Z2,   _}},
-      {{  _,   _,   _,   S,   _,   _,   _,   _,   _,  -S,   _,   _}},
-      {{  _,   _, -Z2,   _,  Z3,   _,   _,   _,  Z2,   _,  Z4,   _}},
-      {{  _,  Y2,   _,   _,   _,  Y3,   _, -Y2,   _,   _,   _,  Y4}},
-      {{ -X,   _,   _,   _,   _,   _,   X,   _,   _,   _,   _,   _}},
-      {{  _, -Y1,   _,   _,   _, -Y2,   _,  Y1,   _,   _,   _, -Y2}},
-      {{  _,   _, -Z1,   _,  Z2,   _,   _,   _,  Z1,   _,  Z2,   _}},
-      {{  _,   _,   _,  -S,   _,   _,   _,   _,   _,   S,   _,   _}},
-      {{  _,   _, -Z2,   _,  Z4,   _,   _,   _,  Z2,   _,  Z3,   _}},
-      {{  _,  Y2,   _,   _,   _,  Y4,   _, -Y2,   _,   _,   _,  Y3}},
-  }};
-  // clang-format on
-
-  return tensor::as_matrix_of_rows(&matrix);
-}
-
-// refer to Cook et. al (2002), "Concepts and applications of Finite Element
-// Analysis", 4th ed., p.26
-template <>
-inline Eigen::Matrix<double, 6, 6, Eigen::RowMajor>
-stiffness_matrix<double, 2>(const Properties<double, 2> &properties,
-                            const double length) {
-  const auto L = length;
-  const auto A = properties.area;
-  const auto E = properties.young_modulus;
-  const auto G = properties.shear_modulus;
-  const auto I_z = properties.area_moment_z;
-  const auto k_y = properties.shear_correction_factor_y;
-
-  const double phi_y = 12 * E * I_z * k_y / A / G / L / L;
-
-  const auto X = A * E / L;
-  const auto Y1 = 12 * E * I_z / (1 + phi_y) / L / L / L;
-  const auto Y2 = 6 * E * I_z / (1 + phi_y) / L / L;
-  const auto Y3 = (4 + phi_y) * E * I_z / (1 + phi_y) / L;
-  const auto Y4 = (2 - phi_y) * E * I_z / (1 + phi_y) / L;
-
-  const auto _ = 0.;
-
-  // clang-format off
-  const tensor::Tensor<double, 6, 6> matrix = {{
-      {{  X,   _,   _,  -X,   _,   _}},
-      {{  _,  Y1,  Y2,   _, -Y1,  Y2}},
-      {{  _,  Y2,  Y3,   _, -Y2,  Y4}},
-      {{ -X,   _,   _,   X,   _,   _}},
-      {{  _, -Y1, -Y2,   _,  Y1, -Y2}},
-      {{  _,  Y2,  Y4,   _, -Y2,  Y3}},
-  }};
-  // clang-format on
-
-  return tensor::as_matrix_of_rows(&matrix);
-}
-
-template <class ValueType_, std::size_t Dimension_>
-Eigen::Matrix<ValueType_, Dimension_ *(Dimension_ + 1),
+template <std::size_t Dimension_>
+Eigen::Matrix<double, Dimension_ *(Dimension_ + 1),
               Dimension_ *(Dimension_ + 1), Eigen::RowMajor>
-rotation_matrix(const tensor::Tensor<ValueType_, Dimension_> &orientation);
-
-// refer to Cook et. al (2002), "Concepts and applications of Finite Element
-// Analysis", 4th ed., p.32
-template <>
-inline Eigen::Matrix<double, 12, 12, Eigen::RowMajor>
-rotation_matrix<double, 3>(const tensor::Tensor<double, 3> &orientation) {
-  // rotation that maps the normalized orientation vector to (1, 0, 0)
-  const auto Lambda = Eigen::Quaternion<double>()
-                          .FromTwoVectors(tensor::as_vector(&orientation),
-                                          Eigen::Vector3d::UnitX())
-                          .normalized()
-                          .toRotationMatrix();
-
-  auto result = Eigen::Matrix<double, 12, 12, Eigen::RowMajor>::Zero().eval();
-
-  result.block(0, 0, 3, 3) = Lambda;
-  result.block(3, 3, 3, 3) = Lambda;
-  result.block(6, 6, 3, 3) = Lambda;
-  result.block(9, 9, 3, 3) = Lambda;
-
-  return result;
-}
-
-// refer to Cook et. al (2002), "Concepts and applications of Finite Element
-// Analysis", 4th ed., p.31
-template <>
-inline Eigen::Matrix<double, 6, 6, Eigen::RowMajor>
-rotation_matrix<double, 2>(const tensor::Tensor<double, 2> &orientation) {
-  const auto normalized =
-      tensor::as_vector(&orientation) / tensor::as_vector(&orientation).norm();
-
-  // rotation that maps the normalized orientation vector to (1, 0)
-  const tensor::Tensor<double, 2, 2> Lambda = {{
-      {{normalized[0], normalized[1]}},
-      {{-normalized[1], normalized[0]}},
-  }};
-
-  auto result = Eigen::Matrix<double, 6, 6, Eigen::RowMajor>::Zero().eval();
-
-  result.block(0, 0, 2, 2) = tensor::as_matrix_of_rows(&Lambda);
-  result(2, 2) = 1.;
-  result.block(3, 3, 2, 2) = tensor::as_matrix_of_rows(&Lambda);
-  result(5, 5) = 1.;
-
-  return result;
-}
+timoshenko_beam_stiffness_matrix(
+    const tensor::Tensor<double, Dimension_> &axis,
+    const Properties<double, Dimension_> &properties) noexcept;
 
 /**
  * @brief Implementation of the closed-form Timoshenko beam element as presented
@@ -213,44 +75,32 @@ rotation_matrix<double, 2>(const tensor::Tensor<double, 2> &orientation) {
  * Analysis", 4th ed., pp.24-32
  */
 template <std::size_t Dimension_>
-struct BeamElement final
-    : ElementBase<BeamElement<Dimension_>, std::size_t, double, 2,
+struct TimoshenkoBeamElement final
+    : ElementBase<TimoshenkoBeamElement<Dimension_>, std::size_t, double, 2,
                   (Dimension_ * (Dimension_ + 1)) / 2> {
 public:
-  using value_type = typename BeamElement::value_type;
-  using size_type = typename BeamElement::size_type;
+  explicit TimoshenkoBeamElement(
+      typename TimoshenkoBeamElement::StiffnessMatrix matrix) noexcept
+      : stiffness_matrix_(std::move(matrix)) {}
 
-  using Vector = tensor::Tensor<value_type, Dimension_>;
-
-  explicit BeamElement(
-      const Vector &element_axis,
-      const Properties<value_type, Dimension_> &properties) noexcept {
-
-    const auto reference_stiffness_matrix =
-        timoshenko::stiffness_matrix<value_type, BeamElement::dimension()>(
-            properties, tensor::as_vector(&element_axis).norm());
-
-    const auto rotation_matrix =
-        timoshenko::rotation_matrix<value_type, BeamElement::dimension()>(
-            element_axis);
-
-    stiffness_matrix_ = rotation_matrix.transpose() *
-                        reference_stiffness_matrix * rotation_matrix;
-  }
-
-  const typename BeamElement::StiffnessMatrix &stiffness_matrix() const {
+  const typename TimoshenkoBeamElement::StiffnessMatrix &
+  stiffness_matrix() const {
     return stiffness_matrix_;
   }
 
-  static constexpr size_type dimension() { return Dimension_; }
+  /**
+   * @brief The dimension of physical space.
+   */
+  static constexpr typename TimoshenkoBeamElement::size_type dimension() {
+    return Dimension_;
+  }
 
 private:
-  typename BeamElement::StiffnessMatrix stiffness_matrix_;
+  typename TimoshenkoBeamElement::StiffnessMatrix stiffness_matrix_;
 };
-} // namespace timoshenko
 
 template <std::size_t Dimension_>
-struct ComputeEnergyTrait<timoshenko::BeamElement<Dimension_>> {
+struct ComputeEnergyTrait<TimoshenkoBeamElement<Dimension_>> {
   template <class Element>
   typename Element::Energy
   operator()(const Element &element,
@@ -263,7 +113,7 @@ struct ComputeEnergyTrait<timoshenko::BeamElement<Dimension_>> {
 };
 
 template <std::size_t Dimension_>
-struct ComputeForcesTrait<timoshenko::BeamElement<Dimension_>> {
+struct ComputeForcesTrait<TimoshenkoBeamElement<Dimension_>> {
   template <class Element>
   typename Element::Forces
   operator()(const Element &element,
@@ -277,7 +127,7 @@ struct ComputeForcesTrait<timoshenko::BeamElement<Dimension_>> {
 };
 
 template <std::size_t Dimension_>
-struct ComputeStiffnessMatrixTrait<timoshenko::BeamElement<Dimension_>> {
+struct ComputeStiffnessMatrixTrait<TimoshenkoBeamElement<Dimension_>> {
   template <class Element>
   typename Element::StiffnessMatrix
   operator()(const Element &element,
diff --git a/elements/test/TimoshenkoBeamElement_Test.cc b/elements/test/TimoshenkoBeamElement_Test.cc
index c90c873e1b73af7ce9f83075cb7fdc221c566411..9f2755ac7b7fb4ad936bf3edc6e07b9be6e9d5fc 100644
--- a/elements/test/TimoshenkoBeamElement_Test.cc
+++ b/elements/test/TimoshenkoBeamElement_Test.cc
@@ -27,43 +27,39 @@ namespace elements {
 namespace {
 
 template <class Element>
-typename Element::Vector create_reference_element_axis() noexcept {
-  auto element_axis = typename Element::Vector();
+tensor::Tensor<double, Element::dimension()>
+create_reference_element_axis() noexcept {
+  auto element_axis = tensor::Tensor<double, Element::dimension()>();
   std::fill(element_axis.begin(), element_axis.end(), 0.);
   element_axis[0] = 1.;
   return element_axis;
 }
 
 template <class Element>
-typename Element::Vector create_rotated_and_stretched_element_axis() noexcept {
-  auto element_axis = typename Element::Vector();
+tensor::Tensor<double, Element::dimension()>
+create_rotated_and_stretched_element_axis() noexcept {
+  auto element_axis = tensor::Tensor<double, Element::dimension()>();
   std::iota(element_axis.begin(), element_axis.end(), 3.);
   return element_axis;
 }
 
-template <class ValueType_, std::size_t Dimension_>
-inline typename timoshenko::Properties<ValueType_, Dimension_>
-create_element_properties() noexcept {
-  exit(1);
-}
+template <std::size_t Dimension_>
+Properties<double, Dimension_> create_element_properties() noexcept;
 
-template <>
-inline typename timoshenko::Properties<double, 3>
-create_element_properties<double, 3>() noexcept {
+template <> Properties<double, 3> create_element_properties<3>() noexcept {
   return {1., 1., 1., 1., 1., 1., 1., 1.};
 }
-template <>
-inline typename timoshenko::Properties<double, 2>
-create_element_properties<double, 2>() noexcept {
+
+template <> Properties<double, 2> create_element_properties<2>() noexcept {
   return {1., 1., 1., 1., 1.};
 }
 
 template <class Element_> struct ReferenceConfiguration {
   using Element = Element_;
   static Element create_element() noexcept {
-    return Element(create_reference_element_axis<Element_>(),
-                   create_element_properties<typename Element_::value_type,
-                                             Element_::dimension()>());
+    return Element(timoshenko_beam_stiffness_matrix(
+        create_reference_element_axis<Element_>(),
+        create_element_properties<Element_::dimension()>()));
   }
 
   static typename Element::Time create_time() noexcept {
@@ -74,9 +70,9 @@ template <class Element_> struct ReferenceConfiguration {
 template <class Element_> struct RotatedAndStretchedConfiguration {
   using Element = Element_;
   static Element create_element() noexcept {
-    return Element(create_rotated_and_stretched_element_axis<Element_>(),
-                   create_element_properties<typename Element_::value_type,
-                                             Element_::dimension()>());
+    return Element(timoshenko_beam_stiffness_matrix(
+        create_rotated_and_stretched_element_axis<Element_>(),
+        create_element_properties<Element_::dimension()>()));
   }
 
   static typename Element::Time create_time() noexcept {
@@ -85,15 +81,15 @@ template <class Element_> struct RotatedAndStretchedConfiguration {
 };
 
 using Configurations =
-    Types<ReferenceConfiguration<timoshenko::BeamElement<2>>,
-          ReferenceConfiguration<timoshenko::BeamElement<3>>,
-          RotatedAndStretchedConfiguration<timoshenko::BeamElement<2>>,
-          RotatedAndStretchedConfiguration<timoshenko::BeamElement<3>>>;
+    Types<ReferenceConfiguration<TimoshenkoBeamElement<2>>,
+          ReferenceConfiguration<TimoshenkoBeamElement<3>>,
+          RotatedAndStretchedConfiguration<TimoshenkoBeamElement<2>>,
+          RotatedAndStretchedConfiguration<TimoshenkoBeamElement<3>>>;
 INSTANTIATE_TYPED_TEST_CASE_P(TimoshenkoBeamElement_Test, Element_Test,
                               Configurations);
 
 struct TimoshenkoBeamElement2D_Test : Test {
-  using Element = timoshenko::BeamElement<2>;
+  using Element = TimoshenkoBeamElement<2>;
   const Element element = ReferenceConfiguration<Element>::create_element();
 };
 
@@ -133,10 +129,11 @@ TEST_F(TimoshenkoBeamElement2D_Test,
 }
 
 struct TimoshenkoBeamElement3D_Test : Test {
-  using Element = timoshenko::BeamElement<3>;
+  using Element = TimoshenkoBeamElement<3>;
   const Element element = ReferenceConfiguration<Element>::create_element();
-  const Element::Vector axis = create_reference_element_axis<Element>();
-  const Element::value_type L = tensor::as_vector(&axis).norm();
+  const tensor::Tensor<double, 3> axis =
+      create_reference_element_axis<Element>();
+  const double L = tensor::as_vector(&axis).norm();
 };
 
 TEST_F(TimoshenkoBeamElement3D_Test, computes_energy_with_axial_displacement) {
@@ -213,12 +210,12 @@ TEST_F(TimoshenkoBeamElement3D_Test,
 }
 
 struct TimoshenkoBeamElement2D_rotated_Test : Test {
-  using Element = timoshenko::BeamElement<2>;
+  using Element = TimoshenkoBeamElement<2>;
   const Element element =
       RotatedAndStretchedConfiguration<Element>::create_element();
-  const Element::Vector axis =
+  const tensor::Tensor<double, 2> axis =
       create_rotated_and_stretched_element_axis<Element>();
-  const Element::value_type L = tensor::as_vector(&axis).norm();
+  const double L = tensor::as_vector(&axis).norm();
 };
 
 TEST_F(TimoshenkoBeamElement2D_rotated_Test,
@@ -260,12 +257,12 @@ TEST_F(TimoshenkoBeamElement2D_rotated_Test,
 }
 
 struct TimoshenkoBeamElement3D_rotated_Test : Test {
-  using Element = timoshenko::BeamElement<3>;
+  using Element = TimoshenkoBeamElement<3>;
   const Element element =
       RotatedAndStretchedConfiguration<Element>::create_element();
-  const Element::Vector axis =
+  const tensor::Tensor<double, 3> axis =
       create_rotated_and_stretched_element_axis<Element>();
-  const Element::value_type L = tensor::as_vector(&axis).norm();
+  const double L = tensor::as_vector(&axis).norm();
 };
 
 TEST_F(TimoshenkoBeamElement3D_rotated_Test,
diff --git a/examples/TimoshenkoBeam.cc b/examples/TimoshenkoBeam.cc
index 2d5f11a50345399ebcde02f4d00e7862c6db29bc..7bebb9467c71a0e7f38367cdba2cf0722ec98cf8 100644
--- a/examples/TimoshenkoBeam.cc
+++ b/examples/TimoshenkoBeam.cc
@@ -78,13 +78,12 @@ constexpr auto vertex_positions = VertexPositions{{
 
 // Now we are ready to select the Timoshenko beam element
 
-using Element = elements::timoshenko::BeamElement<dimension>;
+using Element = elements::TimoshenkoBeamElement<dimension>;
 
 // The Timoshenko beam element comes with a number of geometrical and material
 // related properties. These are stored in the Properties struct
 
-using Properties =
-    elements::timoshenko::Properties<Mesh::value_type, dimension>;
+using Properties = elements::Properties<Mesh::value_type, dimension>;
 
 // Let us define some parameters the linear elastic beam of rectangular cross
 // section
@@ -93,13 +92,11 @@ constexpr Mesh::value_type young_modulus = 1.;
 constexpr Mesh::value_type poisson_ratio = 0.3;
 constexpr Mesh::value_type shear_modulus =
     young_modulus / (2 * (1 + poisson_ratio));
-constexpr Mesh::value_type density = 1.;
 constexpr Mesh::value_type shear_correction_factor_y = 1.2;
 constexpr Mesh::value_type thickness = 1.;
 constexpr Mesh::value_type width = 0.1;
 constexpr Mesh::value_type area = width * thickness;
 constexpr Mesh::value_type area_moment_z = thickness * std::pow(width, 3) / 12.;
-constexpr Mesh::value_type weight = 1.;
 
 // We will assemble e.g. energy using a collection of elements. This is done by
 // the assembler. (The list DefaultFeaturePlugins contain the features (e.g.
@@ -123,9 +120,8 @@ int main(int argc, char **argv) {
         dimension, connectivity, number_of_vertices, dof_per_vertex);
     auto assembler = Assembler();
 
-    Properties properties = {
-        young_modulus, shear_modulus, shear_correction_factor_y, density, area,
-        area_moment_z, weight};
+    Properties properties = {young_modulus, shear_modulus,
+                             shear_correction_factor_y, area, area_moment_z};
 
     // Depending on whether we use MPI, our mesh may be distributed and not all
     // elements are present on this computational node.
@@ -140,7 +136,8 @@ int main(int argc, char **argv) {
           elements::tensor::as_vector(
               &vertex_positions.at(connectivity.at(element.index()).at(0)));
 
-      assembler.emplaceElement(element, element_axis, properties);
+      assembler.emplaceElement(
+          element, timoshenko_beam_stiffness_matrix(element_axis, properties));
     }
 
     // We need to create a solver. We do not use the time, so we can set it to