FESADev/tests/integration/analysis/shell_surface_pressure_nonlinear_test.cpp

#include <gtest/gtest.h>
#include <hdf5.h>

#include <algorithm>
#include <array>
#include <atomic>
#include <chrono>
#include <cmath>
#include <filesystem>
#include <fstream>
#include <iterator>
#include <stdexcept>
#include <string>
#include <vector>

#include "fesa/app/fesa_application.h"
#include "fesa/loads/shell_surface_pressure.h"

namespace {

class TempDirectory {
 public:
  explicit TempDirectory(const std::string& label) {
    static std::atomic<unsigned long long> sequence{0U};
    const auto tick =
        std::chrono::steady_clock::now().time_since_epoch().count();
    path_ =
        std::filesystem::temp_directory_path() /
        ("fesa-surface-pressure-nonlinear-" + label + "-" +
         std::to_string(tick) + "-" + std::to_string(sequence.fetch_add(1U)));
    std::error_code error;
    if (!std::filesystem::create_directory(path_, error) || error) {
      throw std::runtime_error{"Unable to create nonlinear pressure fixture."};
    }
  }

  TempDirectory(const TempDirectory&) = delete;
  TempDirectory& operator=(const TempDirectory&) = delete;

  ~TempDirectory() {
    std::error_code ignored;
    std::filesystem::remove_all(path_, ignored);
  }

  const std::filesystem::path& Path() const noexcept { return path_; }

 private:
  std::filesystem::path path_;
};

class Hdf5File {
 public:
  explicit Hdf5File(const std::filesystem::path& path)
      : handle_{H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)} {
    if (handle_ < 0) {
      throw std::runtime_error{"Unable to open nonlinear pressure results."};
    }
  }

  Hdf5File(const Hdf5File&) = delete;
  Hdf5File& operator=(const Hdf5File&) = delete;

  ~Hdf5File() { (void)H5Fclose(handle_); }

  hid_t Get() const noexcept { return handle_; }

 private:
  hid_t handle_;
};

void WriteText(const std::filesystem::path& path, const std::string& text) {
  std::ofstream stream{path, std::ios::binary | std::ios::trunc};
  stream.write(text.data(), static_cast<std::streamsize>(text.size()));
  if (!stream) {
    throw std::runtime_error{"Unable to write nonlinear pressure input."};
  }
}

std::vector<char> ReadBytes(const std::filesystem::path& path) {
  std::ifstream stream{path, std::ios::binary};
  return {std::istreambuf_iterator<char>{stream},
          std::istreambuf_iterator<char>{}};
}

std::string NonlinearPressureDeck(const bool include_pressure,
                                  const bool collapsed) {
  const std::string surface = include_pressure
                                  ? "*Surface, TYPE=ELEMENT, NAME=Pressure\n"
                                    "Shell-1.1, SPOS\n"
                                  : "";
  const std::string dsload =
      include_pressure ? "*Dsload\nPressure, P, 4.0\n" : "";
  const std::string node_three =
      collapsed ? "3, 1., 0., 0.\n" : "3, 1., 1., 0.\n";
  const std::string constraints =
      "N1, 1, 6\n"
      "N2, 1, 6\n"
      "N3, 2, 6\n"
      "N3, 1, 1, 1.0\n"
      "N4, 1, 6\n";
  return std::string{R"inp(*Part, name=ShellPart
*Node
1, 0., 0., 0.
2, 1., 0., 0.
)inp"} + node_three +
         R"inp(4, 0., 1., 0.
*Element, type=S4
1, 1, 2, 3, 4
*Elset, elset=ShellSet
1
*Shell Section, elset=ShellSet, material=Steel
0.1
*End Part
*Assembly, name=Assembly
*Instance, name=Shell-1, part=ShellPart
*End Instance
)inp" + surface +
         R"inp(*Nset, nset=N1, instance=Shell-1
1
*Nset, nset=N2, instance=Shell-1
2
*Nset, nset=N3, instance=Shell-1
3
*Nset, nset=N4, instance=Shell-1
4
*End Assembly
*Material, name=Steel
*Elastic
1000., 0.25
*Boundary
)inp" + constraints +
         R"inp(*Step, name=Step-1, nlgeom=YES
*Static
0.1, 1., 0.001, 0.1
)inp" + dsload +
         "*End Step\n";
}

std::vector<double> ReadFrameDataset(const Hdf5File& file,
                                     const std::uint32_t frame_index,
                                     const std::string& suffix) {
  const std::string path =
      "/steps/Step-1/frames/" + std::to_string(frame_index) + suffix;
  const hid_t dataset = H5Dopen2(file.Get(), path.c_str(), H5P_DEFAULT);
  if (dataset < 0) {
    throw std::runtime_error{"Unable to open nonlinear pressure dataset."};
  }
  std::vector<double> values(24U);
  const herr_t read = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
                              H5P_DEFAULT, values.data());
  (void)H5Dclose(dataset);
  if (read < 0 ||
      !std::all_of(
          values.begin(), values.end(),
          [](const double value) { return std::isfinite(value); })) {
    throw std::runtime_error{"Unable to read nonlinear pressure dataset."};
  }
  return values;
}

std::array<double, 6> ReadEquilibrium(const Hdf5File& file,
                                      const std::uint32_t frame_index) {
  const std::string path = "/steps/Step-1/frames/" +
                           std::to_string(frame_index) + "/global/equilibrium";
  const hid_t dataset = H5Dopen2(file.Get(), path.c_str(), H5P_DEFAULT);
  if (dataset < 0) {
    throw std::runtime_error{"Unable to open nonlinear equilibrium."};
  }
  std::array<double, 6> values{};
  const herr_t read = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
                              H5P_DEFAULT, values.data());
  (void)H5Dclose(dataset);
  if (read < 0) {
    throw std::runtime_error{"Unable to read nonlinear equilibrium."};
  }
  return values;
}

double SumTranslation(const std::vector<double>& values,
                      const std::size_t component) {
  double result = 0.0;
  for (std::size_t node = 0U; node < 4U; ++node) {
    result += values[node * 6U + component];
  }
  return result;
}

double MaximumAbs(const std::array<double, 6>& values) {
  double result = 0.0;
  for (const double value : values) {
    result = (std::max)(result, std::abs(value));
  }
  return result;
}

fesa::ShellSurfacePressureInput CurrentPressureInput(
    const double node_three_x_displacement) {
  return {
      {{{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {1.0, 1.0, 0.0}, {0.0, 1.0, 0.0}}},
      {{{0.0, 0.0, 0.0},
        {0.0, 0.0, 0.0},
        {node_three_x_displacement, 0.0, 0.0},
        {0.0, 0.0, 0.0}}},
      fesa::ShellSurfaceSide::kSpos,
      4.0,
      1.0,
      fesa::SurfacePressureConfiguration::kCurrent,
      {},
      "unsymmetric"};
}

const Hdf5File RunNonlinearPressure(const std::filesystem::path& input,
                                    const std::filesystem::path& output) {
  testing::internal::CaptureStderr();
  const int exit_code = fesa::FesaApplication{}.Run({input.string()});
  const std::string standard_error = testing::internal::GetCapturedStderr();
  EXPECT_EQ(exit_code, 0) << standard_error;
  EXPECT_TRUE(std::filesystem::exists(output));
  EXPECT_GT(H5Fis_hdf5(output.string().c_str()), 0);
  return Hdf5File{output};
}

}  // namespace

TEST(NonlinearShellSurfacePressure,
     AssemblesFollowerForceFromCurrentTrialGeometry) {
  TempDirectory directory{"current-force"};
  const auto input = directory.Path() / "current-force.inp";
  const auto output = directory.Path() / "current-force_results.h5";
  WriteText(input, NonlinearPressureDeck(true, false));

  const Hdf5File file = RunNonlinearPressure(input, output);
  const auto external =
      ReadFrameDataset(file, 9U, "/nodal/generalized_external_force");

  EXPECT_NEAR(SumTranslation(external, 2U), -6.0, 1.0e-12);
}

TEST(NonlinearShellSurfacePressure,
     FormsResidualAndTangentWithApprovedExternalSigns) {
  TempDirectory directory{"external-sign"};
  const auto input = directory.Path() / "external-sign.inp";
  const auto output = directory.Path() / "external-sign_results.h5";
  WriteText(input, NonlinearPressureDeck(true, false));

  const Hdf5File file = RunNonlinearPressure(input, output);
  const auto external =
      ReadFrameDataset(file, 9U, "/nodal/generalized_external_force");
  const auto residual =
      ReadFrameDataset(file, 9U, "/nodal/generalized_residual");

  EXPECT_NEAR(SumTranslation(external, 2U), -6.0, 1.0e-12);
  EXPECT_NEAR(SumTranslation(residual, 2U), 6.0, 1.0e-12);
}

TEST(NonlinearShellSurfacePressure, AppliesIncrementFactorExactlyOnce) {
  TempDirectory directory{"factor-once"};
  const auto input = directory.Path() / "factor-once.inp";
  const auto output = directory.Path() / "factor-once_results.h5";
  WriteText(input, NonlinearPressureDeck(true, false));

  const Hdf5File file = RunNonlinearPressure(input, output);
  const auto external =
      ReadFrameDataset(file, 0U, "/nodal/generalized_external_force");

  EXPECT_NEAR(SumTranslation(external, 2U), -0.42, 1.0e-12);
}

TEST(NonlinearShellSurfacePressure, UsesFullUnsymmetricTangentSolverPath) {
  TempDirectory directory{"unsymmetric"};
  const auto input = directory.Path() / "unsymmetric.inp";
  const auto output = directory.Path() / "unsymmetric_results.h5";
  WriteText(input, NonlinearPressureDeck(true, false));

  const Hdf5File file = RunNonlinearPressure(input, output);
  const auto external =
      ReadFrameDataset(file, 9U, "/nodal/generalized_external_force");
  const auto residual =
      ReadFrameDataset(file, 9U, "/nodal/generalized_residual");

  EXPECT_NE(SumTranslation(external, 2U), 0.0);
  EXPECT_NE(SumTranslation(residual, 2U), 0.0);

  const auto central =
      fesa::EvaluateShellSurfacePressure(CurrentPressureInput(1.0));
  const auto forward =
      fesa::EvaluateShellSurfacePressure(CurrentPressureInput(1.0 + 1.0e-6));
  const auto backward =
      fesa::EvaluateShellSurfacePressure(CurrentPressureInput(1.0 - 1.0e-6));
  ASSERT_TRUE(central.HasValue());
  ASSERT_TRUE(forward.HasValue());
  ASSERT_TRUE(backward.HasValue());

  bool has_antisymmetric_entry = false;
  for (std::size_t row = 0U; row < 24U; ++row) {
    const double finite_difference =
        (forward.Value().force[row] - backward.Value().force[row]) / 2.0e-6;
    EXPECT_NEAR(finite_difference,
                central.Value().external_tangent[row * 24U + 12U], 1.0e-8)
        << "row=" << row;
    for (std::size_t column = row + 1U; column < 24U; ++column) {
      has_antisymmetric_entry =
          has_antisymmetric_entry ||
          std::abs(central.Value().external_tangent[row * 24U + column] -
                   central.Value().external_tangent[column * 24U + row]) >
              1.0e-12;
    }
  }
  EXPECT_TRUE(has_antisymmetric_entry);
}

TEST(NonlinearShellSurfacePressure,
     ConvergesOpenPatchAndPreservesAntisymmetricEvidence) {
  TempDirectory directory{"open-patch"};
  const auto input = directory.Path() / "open-patch.inp";
  const auto output = directory.Path() / "open-patch_results.h5";
  WriteText(input, NonlinearPressureDeck(true, false));

  const Hdf5File file = RunNonlinearPressure(input, output);
  const auto external =
      ReadFrameDataset(file, 9U, "/nodal/generalized_external_force");

  EXPECT_NEAR(SumTranslation(external, 2U), -6.0, 1.0e-12);
  EXPECT_GT(H5Lexists(file.Get(), "/steps/Step-1/accepted_frames", H5P_DEFAULT),
            0);
}

TEST(NonlinearShellSurfacePressure,
     SatisfiesAcceptedFreeResidualAndFactorOneGlobalBalance) {
  TempDirectory directory{"balance"};
  const auto input = directory.Path() / "balance.inp";
  const auto output = directory.Path() / "balance_results.h5";
  WriteText(input, NonlinearPressureDeck(true, false));

  const Hdf5File file = RunNonlinearPressure(input, output);
  const auto equilibrium = ReadEquilibrium(file, 9U);

  EXPECT_LE(MaximumAbs(equilibrium), 1.0e-8);
}

TEST(NonlinearShellSurfacePressure,
     IsObjectiveUnderRigidMotionAndTracksCurrentArea) {
  fesa::ShellSurfacePressureInput input{
      {{{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {1.0, 1.0, 0.0}, {0.0, 1.0, 0.0}}},
      {{{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, -1.0, 1.0}, {0.0, -1.0, 1.0}}},
      fesa::ShellSurfaceSide::kSpos,
      4.0,
      1.0,
      fesa::SurfacePressureConfiguration::kCurrent,
      {},
      "objectivity"};

  const auto contribution = fesa::EvaluateShellSurfacePressure(input);

  ASSERT_TRUE(contribution.HasValue());
  double force_x = 0.0;
  double force_y = 0.0;
  double force_z = 0.0;
  for (std::size_t node = 0U; node < 4U; ++node) {
    force_x += contribution.Value().force[node * 6U];
    force_y += contribution.Value().force[node * 6U + 1U];
    force_z += contribution.Value().force[node * 6U + 2U];
  }
  EXPECT_NEAR(force_x, 0.0, 1.0e-12);
  EXPECT_NEAR(force_y, 4.0, 1.0e-12);
  EXPECT_NEAR(force_z, 0.0, 1.0e-12);
}

TEST(NonlinearShellSurfacePressure,
     ReducesSerialAndParallelContributionsDeterministically) {
  TempDirectory directory{"deterministic"};
  const auto first_input = directory.Path() / "first.inp";
  const auto first_output = directory.Path() / "first_results.h5";
  const auto second_input = directory.Path() / "second.inp";
  const auto second_output = directory.Path() / "second_results.h5";
  const std::string deck = NonlinearPressureDeck(true, false);
  WriteText(first_input, deck);
  WriteText(second_input, deck);

  const Hdf5File first = RunNonlinearPressure(first_input, first_output);
  const Hdf5File second = RunNonlinearPressure(second_input, second_output);
  const auto first_force =
      ReadFrameDataset(first, 9U, "/nodal/generalized_external_force");
  const auto second_force =
      ReadFrameDataset(second, 9U, "/nodal/generalized_external_force");

  EXPECT_EQ(first_force, second_force);
  EXPECT_NE(SumTranslation(first_force, 2U), 0.0);
}

TEST(NonlinearShellSurfacePressure,
     InvalidTrialRollsBackCommittedFrameAndFinalOutput) {
  TempDirectory directory{"rollback"};
  const auto input = directory.Path() / "rollback.inp";
  const auto output = directory.Path() / "rollback_results.h5";
  WriteText(input, NonlinearPressureDeck(true, true));
  const std::vector<char> previous{'p', 'r', 'i', 'o', 'r'};
  {
    std::ofstream stream{output, std::ios::binary | std::ios::trunc};
    stream.write(previous.data(),
                 static_cast<std::streamsize>(previous.size()));
  }

  const int exit_code = fesa::FesaApplication{}.Run({input.string()});

  EXPECT_NE(exit_code, 0);
  EXPECT_EQ(ReadBytes(output), previous);
}

TEST(NonlinearShellSurfacePressure,
     WritesSchemaOneExternalInternalResidualReactionAndEquilibrium) {
  TempDirectory directory{"schema"};
  const auto input = directory.Path() / "schema.inp";
  const auto output = directory.Path() / "schema_results.h5";
  WriteText(input, NonlinearPressureDeck(true, false));

  const Hdf5File file = RunNonlinearPressure(input, output);
  const std::string root = "/steps/Step-1/frames/9";
  for (const std::string& suffix :
       {"/nodal/generalized_internal_force",
        "/nodal/generalized_external_force", "/nodal/generalized_residual",
        "/nodal/reaction", "/global/equilibrium"}) {
    EXPECT_GT(H5Lexists(file.Get(), (root + suffix).c_str(), H5P_DEFAULT), 0)
        << suffix;
  }
  const auto external =
      ReadFrameDataset(file, 9U, "/nodal/generalized_external_force");
  EXPECT_NEAR(SumTranslation(external, 2U), -6.0, 1.0e-12);
  EXPECT_EQ(H5Lexists(file.Get(), "/surface_pressure", H5P_DEFAULT), 0);
}