compare consistent mass matrix with BE reference

elements/test/TimoshenkoBeamElement_Test.cc

@@ -461,6 +461,78 @@ INSTANTIATE_TYPED_TEST_CASE_P(TimoshenkoBeamElementWithMass_Test,
                               TimoshenkoBeamElementWithMass_Test,
                               MassConfigurations);
 
+template <std::size_t Dimension_>
+TimoshenkoBeamProperties<double, Dimension_>
+create_euler_bernoulli_properties() noexcept {
+  static_assert(Dimension_ == 2 || Dimension_ == 3,
+                "Only dimensions 2 and 3 are supported.");
+  if constexpr (Dimension_ == 2) {
+    return {
+        .young_modulus = 1.,
+        .shear_modulus = 1.,
+        .shear_correction_factor_y = 1.,
+        .area = 1.,
+        .area_moment_z = 0.,
+    };
+  } else {
+    return {
+        .young_modulus = 1.,
+        .shear_modulus = 1.,
+        .shear_correction_factor_y = 1.,
+        .shear_correction_factor_z = 1.,
+        .area = 1.,
+        .area_moment_y = 0.,
+        .area_moment_z = 0.,
+        .polar_moment_x = 0.,
+    };
+  }
+}
+
+struct TimoshenkoBeamElement2D_EulerBernoulliTest : ::testing::Test {
+  using Configuration =
+      ConfigurationWithMass<TimoshenkoBeamElementWithMass<2>,
+                            &create_reference_element_axis<2>,
+                            &create_euler_bernoulli_properties<2>,
+                            &timoshenko_beam_consistent_mass_matrix<2>>;
+  using Element = typename Configuration::Element;
+  Element element = Configuration::create_element();
+
+  const tensor::Tensor<double, 2> axis =
+      create_reference_element_axis<Element::dimension()>();
+  const double L = tensor::as_vector(&axis).norm();
+  const double m = Configuration::calculate_mass();
+};
+
+TEST_F(TimoshenkoBeamElement2D_EulerBernoulliTest,
+       special_case_is_equal_to_euler_bernoulli_mass_matrix) {
+  const auto mass = compute_mass_matrix(element);
+
+  // see Felippa et al (2015), "Mass Matrix Templates: General Description
+  // and 1D Examples", eq. 148, http://dx.doi.org/10.1007/s11831-014-9108-x
+  //
+  // Note that degrees of freedom 0, 3 represent the displacement in axial
+  // direction.
+  const auto reference = [&]() {
+    auto reference = decltype(mass)::Zero().eval();
+    reference(0, 0) = 420. / 3.;
+    reference(0, 3) = reference(3, 0) = 420 / 6.;
+    reference(3, 3) = 420. / 3.;
+    reference(1, 1) = 156.;
+    reference(1, 2) = reference(2, 1) = 22. * L;
+    reference(1, 4) = reference(4, 1) = 54;
+    reference(1, 5) = reference(5, 1) = -13. * L;
+    reference(2, 2) = 4. * L * L;
+    reference(2, 4) = reference(4, 2) = 13. * L;
+    reference(2, 5) = reference(5, 2) = -3. * L * L;
+    reference(4, 4) = 156.;
+    reference(4, 5) = reference(5, 4) = -22. * L;
+    reference(5, 5) = 4 * L * L;
+    return ((m / 420.) * reference).eval();
+  }();
+
+  EXPECT_THAT(mass, IsEigenApprox(reference));
+}
+
 } // namespace
 } // namespace elements
 } // namespace ae108
\ No newline at end of file