#include "fesa/io/hdf5/solid_hdf5_results_writer.h"

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

#include <cfenv>
#include <fstream>
#include <iterator>
#include <limits>
#include <set>

#include "fesa/analysis/linear_frequency_procedure.h"
#include "fesa/analysis/linear_static_procedure.h"
#include "fesa/analysis/newton_iteration_log.h"
#include "fesa/analysis/nonlinear_static_procedure.h"
#include "fesa/analysis/prepared_analysis.h"
#include "fesa/analysis/solid_frequency_procedure.h"
#include "fesa/app/output_paths.h"
#include "fesa/assembly/parallel_for.h"
#include "fesa/io/hdf5/eigen_hdf5_results_writer.h"
#include "fesa/io/hdf5/hdf5_results_writer.h"
#include "fesa/io/hdf5/nonlinear_hdf5_results_writer.h"
#include "fesa/solvers/linear/mkl_pardiso_solver.h"
#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h"
#include "io/hdf5/hdf5_primitives.h"
#include "io/hdf5/solid_hdf5_self_check.h"
#include "solid_output_test_fixture.h"
namespace {
using namespace fesa;
using namespace fesa::hdf5_internal;
struct Actual {
  PreparedAnalysisInput input;
  PreparedAnalysisOutcome outcome;
};
Actual Execute(const solid_output_test::Files& files) {
  auto input = PrepareAnalysisInput(files.input);
  if (!input.HasValue()) throw std::runtime_error("actual input");
  SerialParallelFor parallel;
  MklPardisoSolver linear;
  MklPardisoTangentSolver tangent;
  LinearStaticProcedureExecutor oldlinear(parallel, linear);
  NonlinearStaticProcedureExecutor oldnonlinear(parallel, tangent);
  LinearFrequencyProcedureExecutor oldfrequency;
  SolidStaticProcedureExecutor stat(parallel, linear, tangent);
  SolidFrequencyProcedureExecutor frequency;
  SolidHhtDynamicProcedureExecutor dynamic(parallel, linear, tangent);
  NewtonIterationLogWriter log;
  if (!log.Open(files.log).IsOk()) throw std::runtime_error("log");
  auto result = ExecutePreparedAnalysis(
      input.Value().domain, oldlinear, oldnonlinear, oldfrequency, log, nullptr,
      nullptr, &stat, &frequency, &dynamic);
  if (!log.Close().IsOk() || !result.HasValue())
    throw std::runtime_error("actual execute");
  return {std::move(input.Value()), std::move(result.Value())};
}
SolidStateView View(const Actual& actual) {
  if (const auto* value =
          std::get_if<SolidStaticExecutionResult>(&actual.outcome))
    return std::cref(value->state);
  if (const auto* value = std::get_if<SolidFrequencyState>(&actual.outcome))
    return std::cref(*value);
  return std::cref(std::get<SolidDynamicExecutionResult>(actual.outcome).state);
}
std::vector<char> Bytes(const std::filesystem::path& path) {
  std::ifstream stream(path, std::ios::binary);
  return {std::istreambuf_iterator<char>(stream), {}};
}
void Existing(const std::filesystem::path& path) {
  std::ofstream out(path, std::ios::binary);
  out << "previous final preserved";
}
std::string Attr(hid_t object, const char* name) {
  Hdf5Handle a{H5Aopen(object, name, H5P_DEFAULT), H5Aclose};
  Hdf5Handle type{H5Aget_type(a.Get()), H5Tclose};
  EXPECT_EQ(H5Tget_class(type.Get()), H5T_STRING);
  EXPECT_GT(H5Tis_variable_str(type.Get()), 0);
  EXPECT_EQ(H5Tget_cset(type.Get()), H5T_CSET_UTF8);
  char* pointer = nullptr;
  EXPECT_GE(H5Aread(a.Get(), type.Get(), &pointer), 0);
  std::string value = pointer ? pointer : "";
  H5free_memory(pointer);
  return value;
}
void Compound(hid_t file, const char* path,
              const std::vector<std::string>& names,
              const std::vector<char>& kinds,
              const std::vector<hsize_t>& widths = {}) {
  Hdf5Handle dataset{H5Dopen2(file, path, H5P_DEFAULT), H5Dclose};
  ASSERT_GE(dataset.Get(), 0);
  Hdf5Handle type{H5Dget_type(dataset.Get()), H5Tclose};
  ASSERT_EQ(H5Tget_class(type.Get()), H5T_COMPOUND);
  ASSERT_EQ(H5Tget_nmembers(type.Get()), static_cast<int>(names.size()));
  for (std::size_t i = 0; i < names.size(); ++i) {
    char* name = H5Tget_member_name(type.Get(), static_cast<unsigned>(i));
    ASSERT_NE(name, nullptr);
    EXPECT_EQ(name, names[i]);
    H5free_memory(name);
    Hdf5Handle field{H5Tget_member_type(type.Get(), static_cast<unsigned>(i)),
                     H5Tclose};
    if (!widths.empty() && widths[i] > 1) {
      ASSERT_EQ(H5Tget_class(field.Get()), H5T_ARRAY);
      hsize_t width = 0;
      ASSERT_EQ(H5Tget_array_dims2(field.Get(), &width), 1);
      EXPECT_EQ(width, widths[i]);
      field = Hdf5Handle{H5Tget_super(field.Get()), H5Tclose};
    }
    if (kinds[i] == 's') {
      EXPECT_EQ(H5Tget_class(field.Get()), H5T_STRING);
      EXPECT_GT(H5Tis_variable_str(field.Get()), 0);
      EXPECT_EQ(H5Tget_cset(field.Get()), H5T_CSET_UTF8);
    } else
      EXPECT_GT(H5Tequal(field.Get(), kinds[i] == 'u'   ? H5T_STD_U64LE
                                      : kinds[i] == 'b' ? H5T_STD_U8LE
                                                        : H5T_IEEE_F64LE),
                0);
  }
}

template <class T>
T NumberAttr(hid_t object, const char* name, hid_t file_type,
             hid_t native_type) {
  Hdf5Handle attribute{H5Aopen(object, name, H5P_DEFAULT), H5Aclose};
  Hdf5Handle type{H5Aget_type(attribute.Get()), H5Tclose};
  Hdf5Handle space{H5Aget_space(attribute.Get()), H5Sclose};
  EXPECT_GT(H5Tequal(type.Get(), file_type), 0);
  EXPECT_EQ(H5Sget_simple_extent_type(space.Get()), H5S_SCALAR);
  T value{};
  EXPECT_GE(H5Aread(attribute.Get(), native_type, &value), 0);
  return value;
}
void DatasetContract(hid_t file, const std::string& path,
                     const std::vector<hsize_t>& expected,
                     const char* components, const char* units,
                     const char* location) {
  Hdf5Handle dataset{H5Dopen2(file, path.c_str(), H5P_DEFAULT), H5Dclose};
  ASSERT_GE(dataset.Get(), 0);
  Hdf5Handle type{H5Dget_type(dataset.Get()), H5Tclose};
  EXPECT_GT(H5Tequal(type.Get(), H5T_IEEE_F64LE), 0);
  Hdf5Handle space{H5Dget_space(dataset.Get()), H5Sclose};
  ASSERT_EQ(H5Sget_simple_extent_ndims(space.Get()),
            static_cast<int>(expected.size()));
  std::vector<hsize_t> shape(expected.size());
  H5Sget_simple_extent_dims(space.Get(), shape.data(), nullptr);
  EXPECT_EQ(shape, expected);
  EXPECT_EQ(Attr(dataset.Get(), "component_names"), components);
  EXPECT_EQ(Attr(dataset.Get(), "unit_dimensions"), units);
  EXPECT_EQ(Attr(dataset.Get(), "location"), location);
  EXPECT_EQ(Attr(dataset.Get(), "coordinate_system"), "global-cartesian");
  const auto marker = path.find("/frames/");
  if (marker != std::string::npos) {
    const auto index =
        static_cast<std::uint64_t>(std::stoull(path.substr(marker + 8)));
    Hdf5Handle metadata{H5Gopen2(file, "/metadata", H5P_DEFAULT), H5Gclose};
    const auto procedure = Attr(metadata.Get(), "analysis_procedure");
    const bool dynamic = procedure.find("hht-") == 0;
    const bool nonlinear = procedure == "nonlinear-static-full-newton" ||
                           procedure == "hht-nonlinear-dynamics";
    EXPECT_EQ(Attr(dataset.Get(), "step_name"), "Step-1");
    EXPECT_EQ(NumberAttr<std::uint64_t>(dataset.Get(), "frame_index",
                                        H5T_STD_U64LE, H5T_NATIVE_UINT64),
              index);
    EXPECT_EQ(NumberAttr<std::uint64_t>(dataset.Get(), "increment_index",
                                        H5T_STD_U64LE, H5T_NATIVE_UINT64),
              dynamic                        ? index
              : procedure == "linear-static" ? 1
                                             : index + 1);
    EXPECT_EQ(
        NumberAttr<double>(dataset.Get(), dynamic ? "time" : "load_factor",
                           H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE),
        dynamic                        ? static_cast<double>(index) * .1
        : procedure == "linear-static" ? 1
                                       : static_cast<double>(index + 1) * .5);
    if (path.find("/solid/stress") != std::string::npos)
      EXPECT_EQ(Attr(dataset.Get(), "stress_measure"), "cauchy");
    if (path.find("/solid/strain") != std::string::npos) {
      EXPECT_EQ(Attr(dataset.Get(), "strain_measure"),
                nonlinear ? "green-lagrange" : "infinitesimal");
      EXPECT_EQ(Attr(dataset.Get(), "shear_convention"), "tensor");
    }
    if (path.find("/nodal/acceleration") != std::string::npos)
      EXPECT_EQ(Attr(dataset.Get(), "acceleration_kind"), "physical");
    if (path.find("/nodal/algorithmic_acceleration") != std::string::npos)
      EXPECT_EQ(Attr(dataset.Get(), "acceleration_kind"), "algorithmic");
    if (path.find("/global/strain_energy") != std::string::npos)
      EXPECT_EQ(Attr(dataset.Get(), "definition"),
                nonlinear ? "sum(W*detJ0)" : "0.5*u^T*K0*u");
    if (path.find("/global/kinetic_energy") != std::string::npos)
      EXPECT_EQ(Attr(dataset.Get(), "definition"), "0.5*v^T*M*v");
    if (path.find("/global/external_work") != std::string::npos)
      EXPECT_EQ(Attr(dataset.Get(), "definition"),
                "sum(0.5*(F_old+F_new)*delta_u)");
    if (path.find("/global/energy_balance_error") != std::string::npos)
      EXPECT_EQ(Attr(dataset.Get(), "definition"), "T+U-T0-U0-Wext");
    if (path.find("/global/angular_momentum") != std::string::npos)
      EXPECT_EQ(Attr(dataset.Get(), "origin"), "global0");
  }
}
TEST(SolidHdf5ResultsWriter,
     AllFiveActualApiOutcomesHaveIndependentOfficialSchema) {
  for (unsigned p = 0; p < 5; ++p) {
    SCOPED_TRACE(p);
    solid_output_test::Files files(p);
    auto actual = Execute(files);
    ASSERT_TRUE(WriteSolidFile(files.output, actual.input.domain, View(actual),
                               actual.input.diagnostics)
                    .IsOk());
    Hdf5Handle file{
        H5Fopen(files.output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT),
        H5Fclose};
    ASSERT_GE(file.Get(), 0);
    Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT),
                        H5Gclose};
    EXPECT_EQ(Attr(metadata.Get(), "analysis_procedure"),
              solid_output_test::Procedure(p));
    EXPECT_EQ(Attr(metadata.Get(), "feature_id"), "hex8-structural-analysis");
    EXPECT_EQ(Attr(metadata.Get(), "source_step_name"), "RawSolid");
    EXPECT_EQ(Attr(metadata.Get(), "element_formulation"), "FESA-Hex8");
    EXPECT_EQ(Attr(metadata.Get(), "integration_rule"), "gauss-2x2x2");
    EXPECT_EQ(Attr(metadata.Get(), "mass_formulation"), "consistent-reference");
    EXPECT_EQ(Attr(metadata.Get(), "component_order"), "node-major; UX,UY,UZ");
    Compound(
        file.Get(), "/model/nodes",
        {"internal_node_id", "instance_name", "source_label", "coordinates"},
        {'u', 's', 's', 'd'}, {1, 1, 1, 3});
    Compound(
        file.Get(), "/model/elements",
        {"internal_element_id", "instance_name", "source_label",
         "source_element_type", "internal_formulation", "node_internal_ids",
         "solid_section_internal_id", "material_internal_id"},
        {'u', 's', 's', 's', 's', 'u', 'u', 'u'}, {1, 1, 1, 1, 1, 8, 1, 1});
    Compound(file.Get(), "/model/solid/materials",
             {"internal_material_id", "name", "youngs_modulus", "poisson_ratio",
              "density_present", "density"},
             {'u', 's', 'd', 'd', 'b', 'd'});
    Compound(file.Get(), "/model/solid/sections",
             {"internal_section_id", "source_file", "source_line",
              "source_elset", "material_internal_id"},
             {'u', 's', 'u', 's', 'u'});
    Compound(file.Get(), "/model/solid_integration_points",
             {"internal_element_id", "point_id", "natural_coordinates",
              "reference_coordinates"},
             {'u', 'b', 'd', 'd'}, {1, 1, 3, 3});
    Compound(file.Get(), "/diagnostics",
             {"severity", "code", "file", "line", "keyword", "entity_identity",
              "message"},
             {'s', 's', 's', 'u', 's', 's', 's'});
    Hdf5Handle mask{
        H5Dopen2(file.Get(), "/model/nodal_constraint_mask", H5P_DEFAULT),
        H5Dclose};
    Hdf5Handle masktype{H5Dget_type(mask.Get()), H5Tclose};
    EXPECT_GT(H5Tequal(masktype.Get(), H5T_STD_U8LE), 0);
    EXPECT_EQ(Attr(mask.Get(), "component_names"), "UX,UY,UZ");
    if (p == 2) {
      EXPECT_EQ(H5Lexists(file.Get(), "/steps/Step-1/frames", H5P_DEFAULT), 0);
      Compound(file.Get(), "/steps/Step-1/eigen/modes",
               {"mode_id", "eigenvalue", "classification", "cluster_id"},
               {'u', 'd', 'b', 'u'});
      Compound(file.Get(), "/steps/Step-1/eigen/frequencies",
               {"mode_id", "frequency"}, {'u', 'd'});
      Compound(file.Get(), "/steps/Step-1/eigen/verification/pairs",
               {"sorted_rank", "solver_column", "eigenvalue", "modal_mass",
                "rayleigh_quotient", "eta_r", "eta_norm", "eta_energy"},
               {'u', 'u', 'd', 'd', 'd', 'd', 'd', 'd'});
      DatasetContract(file.Get(), "/steps/Step-1/eigen/vectors", {1, 8, 3},
                      "UX,UY,UZ", "mass-normalized", "nodal");
      for (const auto& quantity :
           std::vector<std::tuple<const char*, const char*, const char*>>{
               {"modes", "eigenvalue", "1/time^2"},
               {"frequencies", "frequency", "1/time"}}) {
        const auto path =
            std::string("/steps/Step-1/eigen/") + std::get<0>(quantity);
        Hdf5Handle dataset{H5Dopen2(file.Get(), path.c_str(), H5P_DEFAULT),
                           H5Dclose};
        EXPECT_EQ(Attr(dataset.Get(), "component_names"),
                  std::get<1>(quantity));
        EXPECT_EQ(Attr(dataset.Get(), "unit_dimensions"),
                  std::get<2>(quantity));
        EXPECT_EQ(Attr(dataset.Get(), "location"), "global");
        EXPECT_EQ(Attr(dataset.Get(), "coordinate_system"), "global-cartesian");
        EXPECT_EQ(H5Aexists(dataset.Get(), "frame_index"), 0);
        EXPECT_EQ(H5Aexists(dataset.Get(), "step_name"), 0);
      }
      const auto& state = std::get<SolidFrequencyState>(actual.outcome);
      ASSERT_EQ(state.eigenpair_checks[0].solver_column, 0u);
    } else {
      const unsigned frames = p == 0 ? 1 : p == 1 ? 2 : 3;
      if (p == 1)
        Compound(file.Get(), "/steps/Step-1/accepted_frames",
                 {"frame_index", "increment_index", "load_factor",
                  "iteration_count", "is_factor_one"},
                 {'u', 'u', 'd', 'u', 'b'});
      if (p >= 3)
        Compound(file.Get(), "/steps/Step-1/accepted_frames",
                 {"frame_index", "increment_index", "time", "iteration_count",
                  "is_initial"},
                 {'u', 'u', 'd', 'u', 'b'});
      if (p == 0)
        EXPECT_EQ(
            H5Lexists(file.Get(), "/steps/Step-1/accepted_frames", H5P_DEFAULT),
            0);
      for (unsigned i = 0; i < frames; ++i) {
        const auto root = "/steps/Step-1/frames/" + std::to_string(i);
        DatasetContract(file.Get(), root + "/nodal/displacement", {8, 3},
                        "U1,U2,U3", "length", "nodal");
        DatasetContract(file.Get(), root + "/nodal/reaction", {8, 3},
                        "RF1,RF2,RF3", "force", "nodal");
        DatasetContract(file.Get(), root + "/nodal/residual", {8, 3},
                        "R1,R2,R3", "force", "nodal");
        DatasetContract(file.Get(), root + "/solid/stress", {1, 8, 6},
                        "S11,S22,S33,S12,S13,S23", "force/length^2",
                        "integration-point");
        DatasetContract(file.Get(), root + "/solid/strain", {1, 8, 6},
                        "E11,E22,E33,E12,E13,E23", "dimensionless",
                        "integration-point");
        DatasetContract(file.Get(), root + "/global/strain_energy", {},
                        "scalar", "force*length", "global");
        if (p >= 3) {
          DatasetContract(file.Get(), root + "/nodal/velocity", {8, 3},
                          "V1,V2,V3", "length/time", "nodal");
          DatasetContract(file.Get(), root + "/nodal/acceleration", {8, 3},
                          "A1,A2,A3", "length/time^2", "nodal");
          DatasetContract(file.Get(), root + "/nodal/algorithmic_acceleration",
                          {8, 3}, "A1,A2,A3", "length/time^2", "nodal");
          for (const auto* energy :
               {"kinetic_energy", "external_work", "energy_balance_error"})
            DatasetContract(file.Get(), root + "/global/" + energy, {},
                            "scalar", "force*length", "global");
          DatasetContract(file.Get(), root + "/global/linear_momentum", {3},
                          "P1,P2,P3", "mass*length/time", "global");
          DatasetContract(file.Get(), root + "/global/angular_momentum", {3},
                          "L1,L2,L3", "mass*length^2/time", "global");
          EXPECT_LE(
              H5Lexists(file.Get(), (root + "/nodal/displacement_low").c_str(),
                        H5P_DEFAULT),
              0);
        }
      }
    }
  }
}
TEST(SolidHdf5ResultsWriter,
     EveryProcedurePreservesPreviousFinalAcrossLifecycleFailures) {
  for (unsigned p = 0; p < 5; ++p) {
    solid_output_test::Files files(p);
    auto actual = Execute(files);
    Existing(files.output);
    const auto previous = Bytes(files.output);
    for (unsigned failure = 0; failure < 5; ++failure) {
      SCOPED_TRACE(p);
      SCOPED_TRACE(failure);
      SolidFileHooks hooks;
      if (failure == 0)
        hooks.before_write = [](hid_t) {
          throw std::runtime_error("injected write");
        };
      if (failure == 1)
        hooks.candidate_close = [](hid_t id) -> herr_t {
          H5Fclose(id);
          return -1;
        };
      if (failure == 2)
        hooks.before_self_check = [](const auto& path) {
          Hdf5Handle file{
              H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT),
              H5Fclose};
          H5Ldelete(file.Get(), "/model/elements", H5P_DEFAULT);
        };
      if (failure == 3)
        hooks.reader_close = [](hid_t id) -> herr_t {
          H5Fclose(id);
          return -1;
        };
      if (failure == 4)
        hooks.before_finalize = [](const auto& path) {
          std::filesystem::remove(path);
        };
      const auto status =
          WriteSolidFile(files.output, actual.input.domain, View(actual),
                         actual.input.diagnostics, hooks);
      EXPECT_FALSE(status.IsOk());
      EXPECT_EQ(Bytes(files.output), previous);
      for (const auto& entry :
           std::filesystem::directory_iterator(files.folder))
        EXPECT_EQ(entry.path().filename().string().find("candidate"),
                  std::string::npos);
    }
  }
}
TEST(SolidHdf5ResultsWriter,
     MalformedWholePairsAndIdentityFailBeforeCandidateCreation) {
  solid_output_test::Files files(3);
  auto actual = Execute(files);
  Existing(files.output);
  const auto previous = Bytes(files.output);
  auto& state = std::get<SolidDynamicExecutionResult>(actual.outcome).state;
  auto& frame = const_cast<SolidDynamicAcceptedFrameState&>(
      state.AcceptedFrames().back());
  const auto original = frame;
  for (unsigned malformed = 0; malformed < 5; ++malformed) {
    SCOPED_TRACE(malformed);
    frame = original;
    if (malformed == 0) frame.low.displacement = Vector{1};
    if (malformed == 1)
      frame.low.displacement[21] = std::numeric_limits<double>::quiet_NaN();
    if (malformed == 2) frame.low.displacement[21] = frame.displacement[21];
    if (malformed == 3) frame.nodal_source_ids[0].source_label_text = "wrong";
    if (malformed == 4) frame.low.residual[21] = 1;
    EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, View(actual),
                                actual.input.diagnostics)
                     .IsOk());
    EXPECT_EQ(Bytes(files.output), previous);
  }
  frame = original;
  const int rounding = std::fegetround();
  ASSERT_EQ(std::fesetround(FE_DOWNWARD), 0);
  const auto status = WriteSolidFile(files.output, actual.input.domain,
                                     View(actual), actual.input.diagnostics);
  std::fesetround(rounding);
  EXPECT_FALSE(status.IsOk());
  EXPECT_EQ(Bytes(files.output), previous);
}
TEST(SolidHdf5ResultsWriter,
     ModalColumnsRemainZeroBasedAndRejectInvalidInventory) {
  solid_output_test::Files files(2);
  auto actual = Execute(files);
  auto& state = std::get<SolidFrequencyState>(actual.outcome);
  ASSERT_EQ(state.eigenpair_checks[0].solver_column, 0u);
  Existing(files.output);
  const auto previous = Bytes(files.output);
  const auto original = state;
  state.eigenpair_checks[0].solver_column = state.eigenpair_checks.size();
  EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, View(actual),
                              actual.input.diagnostics)
                   .IsOk());
  EXPECT_EQ(Bytes(files.output), previous);
  state = original;
  state.backend.column_ids[0] = state.eigenpair_checks.size();
  state.eigenpair_checks[0].solver_column = state.backend.column_ids[0];
  EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, View(actual),
                              actual.input.diagnostics)
                   .IsOk());
}

TEST(SolidHdf5ResultsWriter,
     ReadbackRejectsMalformedTypesShapesValuesAndSourceIds) {
  solid_output_test::Files files(3);
  auto actual = Execute(files);
  Existing(files.output);
  const auto previous = Bytes(files.output);
  for (unsigned malformed = 0; malformed < 9; ++malformed) {
    SCOPED_TRACE(malformed);
    SolidFileHooks hooks;
    hooks.before_self_check = [&](const auto& path) {
      Hdf5Handle file{H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT),
                      H5Fclose};
      const char* result = "/steps/Step-1/frames/2/nodal/displacement";
      if (malformed == 0) H5Ldelete(file.Get(), result, H5P_DEFAULT);
      if (malformed == 1) {
        Hdf5Handle dataset{H5Dopen2(file.Get(), result, H5P_DEFAULT), H5Dclose};
        std::vector<double> values(24, 0);
        values[0] = std::numeric_limits<double>::quiet_NaN();
        H5Dwrite(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
                 H5P_DEFAULT, values.data());
      }
      if (malformed == 2 || malformed == 3) {
        H5Ldelete(file.Get(), result, H5P_DEFAULT);
        hsize_t shape[]{malformed == 2 ? 1u : 8u, 3};
        Hdf5Handle space{H5Screate_simple(2, shape, nullptr), H5Sclose};
        Hdf5Handle data{
            H5Dcreate2(file.Get(), result,
                       malformed == 2 ? H5T_IEEE_F64LE : H5T_IEEE_F32LE,
                       space.Get(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT),
            H5Dclose};
      }
      if (malformed == 4) {
        Hdf5Handle dataset{
            H5Dopen2(file.Get(), "/model/solid_integration_points",
                     H5P_DEFAULT),
            H5Dclose};
        Hdf5Handle memory{H5Tcreate(H5T_COMPOUND, 1), H5Tclose};
        H5Tinsert(memory.Get(), "point_id", 0, H5T_NATIVE_UINT8);
        std::vector<std::uint8_t> values(8, 99);
        H5Dwrite(dataset.Get(), memory.Get(), H5S_ALL, H5S_ALL, H5P_DEFAULT,
                 values.data());
      }
      if (malformed == 5) {
        Hdf5Handle dataset{H5Dopen2(file.Get(), "/model/elements", H5P_DEFAULT),
                           H5Dclose};
        Hdf5Handle memory{H5Tcreate(H5T_COMPOUND, sizeof(std::uint64_t)),
                          H5Tclose};
        H5Tinsert(memory.Get(), "internal_element_id", 0, H5T_NATIVE_UINT64);
        const std::uint64_t value = 99;
        H5Dwrite(dataset.Get(), memory.Get(), H5S_ALL, H5S_ALL, H5P_DEFAULT,
                 &value);
      }
      if (malformed == 6) {
        Hdf5Handle dataset{H5Dopen2(file.Get(), result, H5P_DEFAULT), H5Dclose};
        H5Adelete(dataset.Get(), "unit_dimensions");
        WriteStringAttribute(dataset.Get(), "unit_dimensions", "wrong");
      }
      if (malformed == 7) {
        Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT),
                            H5Gclose};
        H5Adelete(metadata.Get(), "analysis_procedure");
        Hdf5Handle type{H5Tcopy(H5T_C_S1), H5Tclose};
        H5Tset_size(type.Get(), 16);
        Hdf5Handle space{H5Screate(H5S_SCALAR), H5Sclose};
        Hdf5Handle attr{
            H5Acreate2(metadata.Get(), "analysis_procedure", type.Get(),
                       space.Get(), H5P_DEFAULT, H5P_DEFAULT),
            H5Aclose};
      }
      if (malformed == 8) {
        H5Ldelete(file.Get(), "/model/elements", H5P_DEFAULT);
        Hdf5Handle type{H5Tcreate(H5T_COMPOUND, 16), H5Tclose};
        H5Tinsert(type.Get(), "material_internal_id", 0, H5T_STD_U64LE);
        H5Tinsert(type.Get(), "internal_element_id", 8, H5T_STD_U64LE);
        hsize_t one = 1;
        Hdf5Handle space{H5Screate_simple(1, &one, nullptr), H5Sclose};
        Hdf5Handle dataset{
            H5Dcreate2(file.Get(), "/model/elements", type.Get(), space.Get(),
                       H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT),
            H5Dclose};
      }
    };
    EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, View(actual),
                                actual.input.diagnostics, hooks)
                     .IsOk());
    EXPECT_EQ(Bytes(files.output), previous);
  }
}
TEST(SolidHdf5ResultsWriter,
     ReadbackMatchesActualAcceptedValuesAndFullPairObservables) {
  for (unsigned p = 0; p < 5; ++p) {
    SCOPED_TRACE(p);
    solid_output_test::Files files(p);
    auto actual = Execute(files);
    ASSERT_TRUE(WriteSolidFile(files.output, actual.input.domain, View(actual),
                               actual.input.diagnostics)
                    .IsOk());
    Hdf5Handle file{
        H5Fopen(files.output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT),
        H5Fclose};
    const auto compare = [&](const std::string& path,
                             const std::vector<double>& expected) {
      Hdf5Handle dataset{H5Dopen2(file.Get(), path.c_str(), H5P_DEFAULT),
                         H5Dclose};
      ASSERT_GE(dataset.Get(), 0);
      std::vector<double> values(expected.size());
      ASSERT_GE(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
                        H5P_DEFAULT, values.data()),
                0);
      EXPECT_EQ(values, expected);
    };
    const auto vector = [](const Vector& input) {
      std::vector<double> values;
      for (std::size_t i = 0; i < input.Size(); ++i) values.push_back(input[i]);
      return values;
    };
    Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT),
                        H5Gclose};
    EXPECT_EQ(Attr(metadata.Get(), "source_input_identity"),
              SourceInputIdentity(actual.input.domain));
    if (p == 2) {
      const auto& state = std::get<SolidFrequencyState>(actual.outcome);
      std::vector<double> values;
      for (std::size_t i = 0; i < state.vectors.Rows(); ++i)
        for (std::size_t j = 0; j < state.vectors.Columns(); ++j)
          values.push_back(state.vectors(i, j));
      compare("/steps/Step-1/eigen/vectors", values);
    } else {
      const auto frame = [&](const auto& accepted) {
        const auto root =
            "/steps/Step-1/frames/" + std::to_string(accepted.frame_index);
        compare(root + "/nodal/displacement", vector(accepted.displacement));
        compare(root + "/nodal/reaction", vector(accepted.reaction));
        compare(root + "/nodal/residual", vector(accepted.residual));
        std::vector<double> stress, strain;
        for (const auto& point : accepted.solid_results) {
          stress.insert(stress.end(), point.cauchy_stress.begin(),
                        point.cauchy_stress.end());
          strain.insert(strain.end(), point.strain.begin(), point.strain.end());
        }
        compare(root + "/solid/stress", stress);
        compare(root + "/solid/strain", strain);
      };
      if (p < 2)
        for (const auto& accepted :
             std::get<SolidStaticExecutionResult>(actual.outcome)
                 .state.AcceptedFrames())
          frame(accepted);
      else {
        bool nonzero_low = false;
        for (const auto& accepted :
             std::get<SolidDynamicExecutionResult>(actual.outcome)
                 .state.AcceptedFrames()) {
          frame(accepted);
          const auto root =
              "/steps/Step-1/frames/" + std::to_string(accepted.frame_index);
          compare(root + "/nodal/velocity", vector(accepted.velocity));
          compare(root + "/nodal/acceleration", vector(accepted.acceleration));
          compare(root + "/nodal/algorithmic_acceleration",
                  vector(accepted.algorithmic_acceleration));
          compare(root + "/global/strain_energy", {accepted.strain_energy});
          compare(root + "/global/kinetic_energy", {accepted.kinetic_energy});
          compare(root + "/global/external_work", {accepted.external_work});
          compare(root + "/global/energy_balance_error",
                  {accepted.energy_balance_error});
          for (std::size_t i = 0; i < accepted.low.displacement.Size(); ++i)
            nonzero_low = nonzero_low || accepted.low.displacement[i] != 0;
        }
        EXPECT_TRUE(nonzero_low);
      }
    }
  }
}
class FailedTangent final : public TangentSolver {
 public:
  Status FactorizeFresh(const SparseMatrix&) override {
    return Status::Failure(FailureCategory::kSolver,
                           {{Severity::kError,
                             "injected-tangent-failure",
                             {},
                             "",
                             {},
                             "test backend failed"}});
  }
  Status SolveCorrection(const Vector&, Vector&) const override {
    return Status::Failure(FailureCategory::kSolver, {});
  }
};
TEST(SolidPreparedAnalysis,
     FailureRecordsAreLoggedAndLogFailureTakesPrecedence) {
  for (unsigned p : {1u, 3u, 4u}) {
    solid_output_test::Files files(p);
    auto input = PrepareAnalysisInput(files.input);
    ASSERT_TRUE(input.HasValue());
    SerialParallelFor parallel;
    MklPardisoSolver linear;
    FailedTangent tangent;
    LinearStaticProcedureExecutor oldlinear(parallel, linear);
    NonlinearStaticProcedureExecutor oldnonlinear(parallel, tangent);
    LinearFrequencyProcedureExecutor oldfrequency;
    SolidStaticProcedureExecutor stat(parallel, linear, tangent);
    SolidFrequencyProcedureExecutor frequency;
    SolidHhtDynamicProcedureExecutor dynamic(parallel, linear, tangent);
    for (bool closed : {false, true}) {
      SCOPED_TRACE(p);
      SCOPED_TRACE(closed);
      NewtonIterationLogWriter log;
      ASSERT_TRUE(log.Open(files.log).IsOk());
      if (closed) ASSERT_TRUE(log.Close().IsOk());
      const auto result = ExecutePreparedAnalysis(
          input.Value().domain, oldlinear, oldnonlinear, oldfrequency, log,
          nullptr, nullptr, &stat, &frequency, &dynamic);
      ASSERT_FALSE(result.HasValue());
      EXPECT_EQ(result.GetStatus().Category(),
                closed ? FailureCategory::kOutput : FailureCategory::kSolver);
      if (!closed) {
        ASSERT_TRUE(log.Close().IsOk());
        const auto bytes = Bytes(files.log);
        const std::string text(bytes.begin(), bytes.end());
        EXPECT_NE(text.find("injected-tangent-failure"), std::string::npos);
        EXPECT_EQ(text.find("rotation"), std::string::npos);
        EXPECT_NE(text.find("stage="), std::string::npos);
      }
    }
  }
}
TEST(SolidPreparedAnalysis, MissingSolidAdapterFailsWithoutLegacyDispatch) {
  solid_output_test::Files files(0);
  auto input = PrepareAnalysisInput(files.input);
  ASSERT_TRUE(input.HasValue());
  SerialParallelFor parallel;
  MklPardisoSolver linear;
  MklPardisoTangentSolver tangent;
  LinearStaticProcedureExecutor oldlinear(parallel, linear);
  NonlinearStaticProcedureExecutor oldnonlinear(parallel, tangent);
  LinearFrequencyProcedureExecutor oldfrequency;
  NewtonIterationLogWriter log;
  ASSERT_TRUE(log.Open(files.log).IsOk());
  const auto result = ExecutePreparedAnalysis(input.Value().domain, oldlinear,
                                              oldnonlinear, oldfrequency, log);
  ASSERT_FALSE(result.HasValue());
  ASSERT_FALSE(result.GetStatus().Diagnostics().empty());
  EXPECT_EQ(result.GetStatus().Diagnostics()[0].code, "missing-solid-executor");
  EXPECT_TRUE(log.Close().IsOk());
}

TEST(SolidHdf5ResultsWriter, UnsupportedEnvironmentAfterCallbackCannotPublish) {
  solid_output_test::Files files(3);
  auto actual = Execute(files);
  Existing(files.output);
  const auto previous = Bytes(files.output);
  for (unsigned callback = 0; callback < 2; ++callback) {
    for (unsigned control : {0x8000u, 0x40u}) {
      SCOPED_TRACE(callback);
      SCOPED_TRACE(control);
      struct Restore {
        unsigned original = _mm_getcsr();
        ~Restore() { _mm_setcsr(original); }
      } restore;
      SolidFileHooks hooks;
      if (callback == 0)
        hooks.before_write = [&](hid_t) {
          _mm_setcsr(restore.original | control);
        };
      else
        hooks.before_finalize = [&](const auto&) {
          _mm_setcsr(restore.original | control);
        };
      EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain,
                                  View(actual), actual.input.diagnostics, hooks)
                       .IsOk());
      EXPECT_EQ(Bytes(files.output), previous);
    }
  }
}
TEST(SolidPreparedAnalysis, MissingWriterPreservesFinalForEveryTypedOutcome) {
  for (unsigned p = 0; p < 5; ++p) {
    solid_output_test::Files files(p);
    auto actual = Execute(files);
    Existing(files.output);
    const auto previous = Bytes(files.output);
    const auto paths =
        DeriveSingleInputOutputPaths({files.input.string()}, files.folder);
    ASSERT_TRUE(paths.HasValue());
    Hdf5ResultsWriter linear;
    NonlinearHdf5ResultsWriter nonlinear;
    EigenHdf5ResultsWriter eigen;
    const auto status = FinalizePreparedAnalysisResult(
        paths.Value(), actual.input.domain, actual.outcome,
        actual.input.diagnostics, linear, nonlinear, eigen);
    ASSERT_FALSE(status.IsOk());
    ASSERT_FALSE(status.Diagnostics().empty());
    EXPECT_EQ(status.Diagnostics()[0].code, "missing-solid-writer");
    EXPECT_EQ(Bytes(files.output), previous);
  }
}
TEST(SolidNewtonIterationLog,
     InvalidScalarEvidencePermanentlyFailsLogLifecycle) {
  for (unsigned failure = 0; failure < 3; ++failure) {
    solid_output_test::Files files(3);
    NewtonIterationLogWriter log;
    ASSERT_TRUE(log.Open(files.log).IsOk());
    SolidDynamicIterationRecord record;
    record.stage = "failed";
    record.outcome = "failed";
    if (failure == 0)
      record.evidence.correction_metrics[2] =
          std::numeric_limits<double>::infinity();
    if (failure == 1) record.evidence.residual_scale = -1;
    if (failure == 2)
      record.evidence.worst_weighted_terms[0] =
          std::numeric_limits<double>::quiet_NaN();
    EXPECT_FALSE(log.Write(record).IsOk());
    EXPECT_FALSE(log.Close().IsOk());
  }
}
TEST(SolidHdf5ResultsWriter, AbsentStaticDensityIsAnExplicitPlaceholder) {
  solid_output_test::Files files(0);
  auto deck = solid_output_test::Deck(0);
  const auto position = deck.find("*Density\n2\n");
  ASSERT_NE(position, std::string::npos);
  deck.erase(position, std::string("*Density\n2\n").size());
  {
    std::ofstream out(files.input);
    out << deck;
  }
  auto actual = Execute(files);
  ASSERT_TRUE(WriteSolidFile(files.output, actual.input.domain, View(actual),
                             actual.input.diagnostics)
                  .IsOk());
  Hdf5Handle file{
      H5Fopen(files.output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT),
      H5Fclose};
  Hdf5Handle dataset{
      H5Dopen2(file.Get(), "/model/solid/materials", H5P_DEFAULT), H5Dclose};
  struct Density {
    std::uint8_t present;
    double rho;
  } row{};
  Hdf5Handle type{H5Tcreate(H5T_COMPOUND, sizeof(row)), H5Tclose};
  H5Tinsert(type.Get(), "density_present", HOFFSET(Density, present),
            H5T_NATIVE_UINT8);
  H5Tinsert(type.Get(), "density", HOFFSET(Density, rho), H5T_NATIVE_DOUBLE);
  ASSERT_GE(
      H5Dread(dataset.Get(), type.Get(), H5S_ALL, H5S_ALL, H5P_DEFAULT, &row),
      0);
  EXPECT_EQ(row.present, 0);
  EXPECT_EQ(row.rho, 0);
}

TEST(SolidHdf5ResultsWriter,
     ModalPhysicalTableAttributeCorruptionPreservesFinal) {
  solid_output_test::Files files(2);
  auto actual = Execute(files);
  Existing(files.output);
  const auto previous = Bytes(files.output);
  for (const auto* quantity : {"modes", "frequencies"}) {
    SolidFileHooks hooks;
    hooks.before_self_check = [&](const auto& path) {
      Hdf5Handle file{H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT),
                      H5Fclose};
      const auto dataset_path = std::string("/steps/Step-1/eigen/") + quantity;
      Hdf5Handle dataset{
          H5Dopen2(file.Get(), dataset_path.c_str(), H5P_DEFAULT), H5Dclose};
      H5Adelete(dataset.Get(), "unit_dimensions");
      WriteStringAttribute(dataset.Get(), "unit_dimensions", "wrong");
    };
    EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, View(actual),
                                actual.input.diagnostics, hooks)
                     .IsOk());
    EXPECT_EQ(Bytes(files.output), previous);
  }
}

TEST(SolidHdf5ResultsWriter,
     UnselectedBackendPairsCannotPublishInvalidSpectrumOrDimension) {
  solid_output_test::Files files(2);
  auto deck = solid_output_test::Deck(2);
  const auto position = deck.find("Tip,2,3\n");
  ASSERT_NE(position, std::string::npos);
  deck.erase(position, std::string("Tip,2,3\n").size());
  {
    std::ofstream out(files.input);
    out << deck;
  }
  auto actual = Execute(files);
  auto& state = std::get<SolidFrequencyState>(actual.outcome);
  ASSERT_EQ(state.modes.size(), 1u);
  ASSERT_EQ(state.free_dof_count, 3u);
  ASSERT_GT(state.eigenpair_checks.size(), state.modes.size());
  ASSERT_LE(state.eigenpair_checks.size(), state.free_dof_count);
  const auto original = state;
  Existing(files.output);
  const auto previous = Bytes(files.output);
  for (unsigned failure = 0; failure < 3; ++failure) {
    SCOPED_TRACE(failure);
    state = original;
    if (failure == 0)
      state.eigenpair_checks[1].eigenvalue = -1e-6 * state.classification_scale;
    if (failure == 1)
      state.eigenpair_checks[1].eigenvalue =
          state.eigenpair_checks[0].eigenvalue * .5;
    if (failure < 2)
      state.backend.eigenvalues[1] = state.eigenpair_checks[1].eigenvalue;
    if (failure == 2) {
      while (state.eigenpair_checks.size() <= state.free_dof_count) {
        auto pair = state.eigenpair_checks.back();
        pair.sorted_rank = state.eigenpair_checks.size() + 1;
        pair.solver_column = state.eigenpair_checks.size();
        state.eigenpair_checks.push_back(pair);
        state.backend.eigenvalues.push_back(pair.eigenvalue);
        state.backend.column_ids.push_back(pair.solver_column);
      }
    }
    const auto status = WriteSolidFile(files.output, actual.input.domain,
                                       View(actual), actual.input.diagnostics);
    ASSERT_FALSE(status.IsOk());
    ASSERT_FALSE(status.Diagnostics().empty());
    const auto& message = status.Diagnostics()[0].message;
    EXPECT_NE(message.find(failure == 0   ? "Significant negative"
                           : failure == 1 ? "not sorted"
                                          : "inventory"),
              std::string::npos);
    EXPECT_EQ(Bytes(files.output), previous);
  }
}
}  // namespace