use Tensor to define 3D stiffness matrix

elements/src/include/ae108/elements/TimoshenkoBeamElement.h

@@ -79,7 +79,6 @@ 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;
@@ -90,45 +89,39 @@ stiffness_matrix<double, 3>(const Properties<double, 3> &properties,
   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;
-
-  auto K = Eigen::Matrix<double, 12, 12, Eigen::RowMajor>::Zero().eval();
+  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 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;
 
-  K(0, 0) = K(6, 6) = X;
-  K(1, 1) = K(7, 7) = Y1;
-  K(5, 5) = K(11, 11) = Y3;
-  K(2, 2) = K(8, 8) = Z1;
-  K(3, 3) = K(9, 9) = S;
-  K(4, 4) = K(10, 10) = Z3;
-
-  K(6, 0) = -X;
-  K(7, 1) = -Y1;
-  K(5, 1) = Y2;
-  K(11, 1) = Y2;
-  K(7, 5) = -Y2;
-  K(11, 7) = -Y2;
-  K(11, 5) = Y4;
-  K(8, 2) = -Z1;
-  K(4, 2) = -Z2;
-  K(10, 2) = -Z2;
-  K(8, 4) = Z2;
-  K(10, 8) = Z2;
-  K(10, 4) = Z4;
-  K(9, 3) = -S;
-
-  return properties.weight * K.template selfadjointView<Eigen::Upper>();
+  const auto _ = 0.;
+
+  const tensor::Tensor<double, 12, 12> matrix = {{
+      {{X, 0., 0., 0., 0., 0., -X, 0., 0., 0., 0., 0.}},
+      {{_, Y1, 0., 0., 0., Y2, 0., -Y1, 0., 0., 0., Y2}},
+      {{_, _, Z1, 0., -Z2, 0., 0., 0., -Z1, 0., -Z2, 0.}},
+      {{_, _, _, S, 0., 0., 0., 0., 0., -S, 0., 0.}},
+      {{_, _, _, _, Z3, 0., 0., 0., Z2, 0., Z4, 0.}},
+      {{_, _, _, _, _, Y3, 0. - Y2, 0., 0., 0., Y4}},
+      {{_, _, _, _, _, _, X, 0., 0., 0., 0., 0.}},
+      {{_, _, _, _, _, _, _, Y1, 0., 0., 0., -Y2}},
+      {{_, _, _, _, _, _, _, _, Z1, 0., Z2, 0.}},
+      {{_, _, _, _, _, _, _, _, _, S, 0., 0.}},
+      {{_, _, _, _, _, _, _, _, _, _, Z3, 0.}},
+      {{_, _, _, _, _, _, _, _, _, _, _, Y3}},
+  }};
+
+  return properties.weight * tensor::as_matrix_of_rows(&matrix)
+                                 .template selfadjointView<Eigen::Upper>();
 }
 
 // refer to Cook et. al (2002), "Concepts and applications of Finite Element
@@ -137,7 +130,6 @@ 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;