FESADev/tests/integration/analysis/nonlinear_static_procedure_test.cpp

#include "fesa/app/fesa_application.h"

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

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

#include "fesa/elements/finite_rotation.h"
#include "fesa/analysis/prepared_analysis.h"
#include "fesa/analysis/nonlinear_static_procedure.h"
#include "fesa/analysis/newton_iteration_log.h"
#include "fesa/assembly/parallel_for.h"
#include "fesa/elements/element_factory.h"
#include "fesa/elements/mitc4_shell.h"
#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h"
#include "fesa/math/vector3.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-step10-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 procedure 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_;
};

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 procedure input."};
  }
}

std::string AllConstrainedNonlinearShellDeck() {
  return R"inp(*Part, name=ShellPart
*Node
1, 0., 0., 0.
2, 1., 0., 0.
3, 1., 1., 0.
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
*Nset, nset=All, instance=Shell-1
1, 2, 3, 4
*End Assembly
*Material, name=Steel
*Elastic
1000., 0.25
*Boundary
All, 1, 6
*Step, name=Step-1, nlgeom=YES
*Static
0.1, 1., 0.001, 0.1
*End Step
)inp";
}

std::string NonzeroPrescribedNonlinearShellDeck() {
  return R"inp(*Part, name=ShellPart
*Node
1, 0., 0., 0.
2, 1., 0., 0.
3, 1., 1., 0.
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
*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
N1, 1, 6
N2, 1, 1, 0.2
N2, 2, 3
N2, 4, 4, 0.3
N2, 5, 6
N3, 1, 6
N4, 1, 6
*Step, name=Step-1, nlgeom=YES
*Static
0.1, 1., 0.001, 0.1
*Cload
N2, 3, -3.
N2, 4, 2.
*End Step
)inp";
}

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 results HDF5."};
    }
  }

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

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

  hid_t Get() const noexcept { return handle_; }

 private:
  hid_t handle_;
};

std::string ReadStringAttribute(const hid_t object, const char* name) {
  const hid_t attribute = H5Aopen(object, name, H5P_DEFAULT);
  const hid_t type = attribute < 0 ? -1 : H5Aget_type(attribute);
  if (attribute < 0 || type < 0) {
    if (type >= 0) {
      (void)H5Tclose(type);
    }
    if (attribute >= 0) {
      (void)H5Aclose(attribute);
    }
    throw std::runtime_error{"Unable to read nonlinear metadata attribute."};
  }
  char* value = nullptr;
  const herr_t read = H5Aread(attribute, type, &value);
  const std::string result = read >= 0 && value != nullptr ? value : "";
  if (value != nullptr) {
    (void)H5free_memory(value);
  }
  (void)H5Tclose(type);
  (void)H5Aclose(attribute);
  if (read < 0) {
    throw std::runtime_error{"Unable to read nonlinear metadata text."};
  }
  return result;
}

std::vector<double> ReadNodalDataset(const Hdf5File& file,
                                     const std::string& path) {
  const hid_t dataset = H5Dopen2(file.Get(), path.c_str(), H5P_DEFAULT);
  if (dataset < 0) {
    throw std::runtime_error{"Unable to open nonlinear nodal 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 nodal dataset."};
  }
  return values;
}

std::vector<std::uint8_t> ReadConstraintMask(const Hdf5File& file) {
  const hid_t dataset =
      H5Dopen2(file.Get(), "/model/nodal_constraint_mask", H5P_DEFAULT);
  if (dataset < 0) {
    throw std::runtime_error{"Unable to open nonlinear constraint mask."};
  }
  std::vector<std::uint8_t> values(24U);
  const herr_t read = H5Dread(dataset, H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL,
                              H5P_DEFAULT, values.data());
  (void)H5Dclose(dataset);
  if (read < 0 || !std::all_of(values.begin(), values.end(),
                                [](const std::uint8_t value) {
                                  return value <= 1U;
                                })) {
    throw std::runtime_error{"Unable to read nonlinear constraint mask."};
  }
  return values;
}

fesa::Vector3 SpatialMoment(const fesa::Vector3& phi,
                            const fesa::Vector3& generalized) {
  // Independent closed-form action, used only as a test oracle.
  const double angle = phi.Norm();
  const double square = angle * angle;
  const double coefficient = angle < 1e-4 ?
      1.0 / 12.0 + square / 720.0 + square * square / 30240.0 +
          square * square * square / 1209600.0 :
      (1.0 - 0.5 * angle / std::tan(0.5 * angle)) / square;
  return generalized + 0.5 * phi.Cross(generalized) +
         coefficient * phi.Cross(phi.Cross(generalized));
}

std::array<double, 6> ReconstructGlobalBalance(const Hdf5File& file) {
  const std::string root = "/steps/Step-1/frames/9/nodal/";
  const auto displacement = ReadNodalDataset(file, root + "displacement");
  const auto external =
      ReadNodalDataset(file, root + "generalized_external_force");
  const auto residual =
      ReadNodalDataset(file, root + "generalized_residual");
  const auto mask = ReadConstraintMask(file);
  const std::array<fesa::Vector3, 4> initial_positions{
      fesa::Vector3{0.0, 0.0, 0.0}, fesa::Vector3{1.0, 0.0, 0.0},
      fesa::Vector3{1.0, 1.0, 0.0}, fesa::Vector3{0.0, 1.0, 0.0}};
  fesa::Vector3 applied_force_resultant;
  fesa::Vector3 applied_moment_resultant;
  fesa::Vector3 reaction_force_resultant;
  fesa::Vector3 reaction_moment_resultant;
  for (std::size_t node = 0U; node < initial_positions.size(); ++node) {
    const std::size_t offset = node * 6U;
    const fesa::Vector3 current_position =
        initial_positions[node] +
        fesa::Vector3{displacement[offset], displacement[offset + 1U],
                      displacement[offset + 2U]};
    const fesa::Vector3 applied_force{external[offset], external[offset + 1U],
                                      external[offset + 2U]};
    const fesa::Vector3 reaction_force{
        mask[offset] == 1U ? residual[offset] : 0.0,
        mask[offset + 1U] == 1U ? residual[offset + 1U] : 0.0,
        mask[offset + 2U] == 1U ? residual[offset + 2U] : 0.0};
    const fesa::Vector3 phi{displacement[offset + 3U],
                             displacement[offset + 4U],
                             displacement[offset + 5U]};
    const fesa::Vector3 applied_moment = SpatialMoment(
        phi, {external[offset + 3U], external[offset + 4U],
              external[offset + 5U]});
    const fesa::Vector3 reaction_moment = SpatialMoment(
        phi, {mask[offset + 3U] == 1U ? residual[offset + 3U] : 0.0,
              mask[offset + 4U] == 1U ? residual[offset + 4U] : 0.0,
              mask[offset + 5U] == 1U ? residual[offset + 5U] : 0.0});
    applied_force_resultant = applied_force_resultant + applied_force;
    applied_moment_resultant =
        applied_moment_resultant + current_position.Cross(applied_force) +
        applied_moment;
    reaction_force_resultant = reaction_force_resultant + reaction_force;
    reaction_moment_resultant =
        reaction_moment_resultant + current_position.Cross(reaction_force) +
        reaction_moment;
  }
  const fesa::Vector3 force_balance =
      applied_force_resultant + reaction_force_resultant;
  const fesa::Vector3 moment_balance =
      applied_moment_resultant + reaction_moment_resultant;
  return {force_balance.X(), force_balance.Y(), force_balance.Z(),
          moment_balance.X(), moment_balance.Y(), moment_balance.Z()};
}

// NSNR-PROC-004, NSNR-PROC-007, and NSNR-PROC-010.
TEST(NonlinearStaticProcedure,
     AllConstrainedNlgeomRunWritesSchemaV1AfterAcceptedFrames) {
  TempDirectory directory{"schema-v1"};
  const auto input = directory.Path() / "shellnl.inp";
  const auto output = directory.Path() / "shellnl_results.h5";
  WriteText(input, AllConstrainedNonlinearShellDeck());

  testing::internal::CaptureStderr();
  const int exit_code = fesa::FesaApplication{}.Run({input.string()});
  const std::string standard_error = testing::internal::GetCapturedStderr();

  ASSERT_EQ(exit_code, 0) << standard_error;
  ASSERT_TRUE(std::filesystem::exists(output));
  ASSERT_GT(H5Fis_hdf5(output.string().c_str()), 0);
  const Hdf5File file{output};
  const hid_t metadata = H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT);
  ASSERT_GE(metadata, 0);
  EXPECT_EQ(ReadStringAttribute(metadata, "feature_id"),
            "nonlinear-static-newton-raphson");
  (void)H5Gclose(metadata);
  EXPECT_GT(H5Lexists(file.Get(), "/steps/Step-1/accepted_frames", H5P_DEFAULT),
            0);
}

// NSNR-PROC-019.
TEST(NonlinearStaticProcedure,
     NonzeroPrescribedMotionWritesCurrentArmSpatialAppliedPlusReactionBalance) {
  TempDirectory directory{"current-arm-balance"};
  const auto input = directory.Path() / "shellnl-current-arm.inp";
  const auto output = directory.Path() / "shellnl-current-arm_results.h5";
  WriteText(input, NonzeroPrescribedNonlinearShellDeck());

  testing::internal::CaptureStderr();
  const int exit_code = fesa::FesaApplication{}.Run({input.string()});
  const std::string standard_error = testing::internal::GetCapturedStderr();

  ASSERT_EQ(exit_code, 0) << standard_error;
  const Hdf5File file{output};
  const hid_t dataset = H5Dopen2(
      file.Get(), "/steps/Step-1/frames/9/global/equilibrium", H5P_DEFAULT);
  ASSERT_GE(dataset, 0);
  std::array<double, 6> actual{};
  ASSERT_GE(H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
                    H5P_DEFAULT, actual.data()),
            0);
  (void)H5Dclose(dataset);
  const auto expected = ReconstructGlobalBalance(file);
  for (std::size_t component = 0U; component < actual.size(); ++component) {
    EXPECT_NEAR(actual[component], expected[component], 1.0e-10)
        << "component=" << component;
  }
}

}  // namespace

TEST(Mitc4DrillingEvidence, NonlinearAndPressureFramesRetainSignedTorqueOracle) {
  for (const bool pressure : {false, true}) {
    SCOPED_TRACE(pressure);
    TempDirectory directory{"drilling"};
    auto deck = AllConstrainedNonlinearShellDeck();
    const auto boundary = deck.find("All, 1, 6");
    ASSERT_NE(boundary, std::string::npos);
    deck.replace(boundary, std::string{"All, 1, 6"}.size(), "All, 1, 5\nAll, 6, 6, 0.125");
    if (pressure) {
      const auto assembly_end = deck.find("*End Assembly");
      deck.insert(assembly_end, "*Surface, TYPE=ELEMENT, NAME=Pressure\nShell-1.1, SPOS\n");
      deck.insert(deck.find("*End Step"), "*Dsload\nPressure, P, 4.0\n");
    }
    const auto path = directory.Path() / "drilling.inp";
    WriteText(path, deck);
    auto prepared = fesa::PrepareAnalysisInput(path);
    ASSERT_TRUE(prepared.HasValue());
    const auto& domain = prepared.Value().domain;
    auto element = fesa::ElementFactory{}.Create(domain.Elements()[0U], domain);
    ASSERT_TRUE(element.HasValue());
    const auto* shell = dynamic_cast<const fesa::Mitc4Shell*>(element.Value().get());
    ASSERT_NE(shell, nullptr);
    const auto stiffness = shell->Stiffness();
    ASSERT_TRUE(stiffness.HasValue());
    double trace = 0.0;
    for (std::size_t node = 0U; node < 4U; ++node)
      for (std::size_t c = 3U; c < 5U; ++c)
        trace += stiffness.Value().physical_local20(5U * node + c, 5U * node + c);
    const double kd = 1e-3 * trace / 8.0;
    fesa::SerialParallelFor serial;
    fesa::MklPardisoTangentSolver tangent;
    fesa::NewtonIterationLogWriter log;
    ASSERT_TRUE(log.Open(directory.Path() / "drilling.log").IsOk());
    auto result = fesa::NonlinearStaticProcedureExecutor{serial, tangent}.Execute(domain, log);
    ASSERT_TRUE(log.Close().IsOk());
    ASSERT_TRUE(result.HasValue()) << (result.GetStatus().Diagnostics().empty() ? "" : result.GetStatus().Diagnostics()[0U].message);
    ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 10U);
    for (const auto& frame : result.Value().state.AcceptedFrames()) {
      ASSERT_TRUE(frame.mitc4_drilling);
      EXPECT_EQ(frame.mitc4_drilling->mode, fesa::Mitc4DrillingMode::kDirectSpring);
      const double torque = 4.0 * kd * 0.125 * frame.load_factor;
      EXPECT_NEAR(frame.global_equilibrium[5U], torque, 1e-13);
      EXPECT_NEAR(frame.mitc4_drilling->drilling_moment[2U], torque, 1e-13);
      for (const double value : frame.mitc4_drilling->equilibrium_with_drilling)
        EXPECT_NEAR(value, 0.0, 1e-13);
      EXPECT_NEAR(frame.physical_strain_energy, 0.0, 1e-20);
      if (pressure) EXPECT_NE(frame.external_force[2U], 0.0);
    }
    ASSERT_TRUE(result.Value().state.CommittedState().mitc4_drilling);
    EXPECT_EQ(result.Value().state.CommittedState().mitc4_drilling->drilling_moment,
              result.Value().state.AcceptedFrames().back().mitc4_drilling->drilling_moment);
  }
}