FESADev/tests/unit/analysis/eigen_mass_assembly_test.cpp

#include "fesa/analysis/linear_frequency_procedure.h"

#include <gtest/gtest.h>

#include <algorithm>
#include <array>
#include <atomic>
#include <chrono>
#include <cmath>
#include <filesystem>
#include <fstream>
#include <sstream>
#include <string>

#include "fesa/analysis/prepared_analysis.h"
#include "fesa/elements/element_factory.h"
#include "fesa/elements/mitc4_shell.h"
#include "fesa/model/domain.h"

namespace {

std::string StatusText(const fesa::Status& status) {
  std::string result;
  for (const auto& diagnostic : status.Diagnostics())
    result += diagnostic.code + " " + diagnostic.message + "\n";
  return result;
}

fesa::Result<fesa::Domain> SharedNodeDomain(const bool reverse_source_order) {
  std::ostringstream deck;
  deck << R"inp(*Part,name=Panel
*Node
1,0,0,0
2,1,0,0
3,2,0,0
4,0,1,0
5,1,1,0
6,2,1,0
*Element,type=S4
)inp";
  if (reverse_source_order)
    deck << "20,2,3,6,5\n10,1,2,5,4\n";
  else
    deck << "10,1,2,5,4\n20,2,3,6,5\n";
  deck << R"inp(*Elset,elset=Thin
10
*Elset,elset=Thick
20
*Shell Section,elset=Thin,material=Light
0.5
*Shell Section,elset=Thick,material=Heavy
1.0
*End Part
*Assembly,name=Assembly
*Instance,name=Panel-1,part=Panel
*End Instance
*Nset,nset=Fixed,instance=Panel-1
1,2,3,4,6
*Nset,nset=Shared,instance=Panel-1
5
*End Assembly
*Material,name=Light
*Elastic
12,0
*Density
2
*Material,name=Heavy
*Elastic
12,0
*Density
5
*Boundary
Fixed,1,6
Shared,1,3
Shared,5,6
*Step,name=Step-2,nlgeom=NO,perturbation
*Frequency
1
*End Step
)inp";
  static std::atomic<unsigned long long> sequence{0U};
  const auto path = std::filesystem::temp_directory_path() /
                    ("fesa-eigen-mass-" +
                     std::to_string(std::chrono::steady_clock::now()
                                        .time_since_epoch()
                                        .count()) +
                     "-" + std::to_string(sequence.fetch_add(1U)) + ".inp");
  {
    std::ofstream stream{path};
    stream << deck.str();
  }
  auto prepared = fesa::PrepareAnalysisInput(path);
  std::error_code ignored;
  std::filesystem::remove(path, ignored);
  if (!prepared.HasValue())
    return fesa::Result<fesa::Domain>::Failure(prepared.GetStatus());
  return fesa::Result<fesa::Domain>::Success(
      std::move(prepared.Value().domain));
}

struct SharedDofPencil {
  double stiffness{0.0};
  double mass{0.0};
};

fesa::Result<SharedDofPencil> DirectSharedDofPencil(
    const fesa::Domain& domain) {
  const fesa::ElementFactory factory;
  SharedDofPencil result;
  for (std::size_t index = 0U; index < domain.ShellElements().Size(); ++index) {
    const auto& definition = domain.ShellElements()[index];
    const auto runtime = factory.Create(definition, domain);
    if (!runtime.HasValue())
      return fesa::Result<SharedDofPencil>::Failure(runtime.GetStatus());
    const auto* shell =
        dynamic_cast<const fesa::Mitc4Shell*>(runtime.Value().get());
    if (shell == nullptr)
      return fesa::Result<SharedDofPencil>::Failure(fesa::Status::Failure(
          fesa::FailureCategory::kModel,
          {{fesa::Severity::kError, "invalid-test-shell", {}, "ELEMENT", {},
            "Frequency mass fixture did not create an MITC4 shell."}}));
    const auto local = std::find(definition.NodeIndices().begin(),
                                 definition.NodeIndices().end(), 4U);
    if (local == definition.NodeIndices().end())
      return fesa::Result<SharedDofPencil>::Failure(fesa::Status::Failure(
          fesa::FailureCategory::kModel,
          {{fesa::Severity::kError, "invalid-test-shell", {}, "ELEMENT", {},
            "Frequency mass fixture lost its shared node."}}));
    const auto mass =
        shell->Mass(definition.MaterialIndex() == 0U ? 2.0 : 5.0);
    const auto stiffness = shell->Stiffness();
    if (!mass.HasValue())
      return fesa::Result<SharedDofPencil>::Failure(mass.GetStatus());
    if (!stiffness.HasValue())
      return fesa::Result<SharedDofPencil>::Failure(stiffness.GetStatus());
    const std::size_t dof =
        6U * static_cast<std::size_t>(
                 std::distance(definition.NodeIndices().begin(), local)) +
        3U;
    result.mass += mass.Value().lumped_global24(dof, dof);
    result.stiffness += stiffness.Value().stabilized_global24(dof, dof);
  }
  return fesa::Result<SharedDofPencil>::Success(result);
}

// E2-C06: sections and densities differ across the shared source node.  The
// one free physical rotation makes the executor's reduced pencil exactly the
// independently scattered local diagonal sum.
TEST(EigenMassAssembly, ScattersHeterogeneousSharedNodeMassInSourceOrder) {
  constexpr double analytic_shared_rotary_mass = 7.0 / 64.0;
  std::array<double, 2> eigenvalues{};
  for (std::size_t order = 0U; order < eigenvalues.size(); ++order) {
    const auto domain = SharedNodeDomain(order != 0U);
    ASSERT_TRUE(domain.HasValue()) << StatusText(domain.GetStatus());
    const auto expected = DirectSharedDofPencil(domain.Value());
    ASSERT_TRUE(expected.HasValue()) << StatusText(expected.GetStatus());
    EXPECT_LE(std::abs(expected.Value().mass - analytic_shared_rotary_mass) /
                  analytic_shared_rotary_mass,
              1e-12);
    const auto state =
        fesa::LinearFrequencyProcedureExecutor{}.Execute(domain.Value());
    ASSERT_TRUE(state.HasValue()) << StatusText(state.GetStatus());
    ASSERT_EQ(state.Value().free_dof_count, 1U);
    ASSERT_EQ(state.Value().modes.size(), 1U);
    const double expected_eigenvalue =
        expected.Value().stiffness / analytic_shared_rotary_mass;
    EXPECT_LE(std::abs(state.Value().modes[0U].eigenvalue - expected_eigenvalue) /
                  (std::abs(state.Value().modes[0U].eigenvalue) +
                   std::abs(expected_eigenvalue)),
              1e-12);
    const double amplitude = state.Value().vectors(0U, 6U * 4U + 3U);
    EXPECT_LE(std::abs(amplitude * amplitude * analytic_shared_rotary_mass -
                       1.0),
              1e-12);
    eigenvalues[order] = state.Value().modes[0U].eigenvalue;
  }
  EXPECT_NEAR(eigenvalues[0U], eigenvalues[1U], 1e-12);
}

}  // namespace