#include "io/hdf5/hdf5_result_writer.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <filesystem>
#include <limits>
#include <set>
#include <string>
#include <utility>
#include <vector>

#include "fesa/analysis/analysis_model.h"
#include "fesa/fem/dof_manager.h"
#include "fesa/math/vector3.h"

namespace fesa::hdf5_internal {

struct StressWriteRow {
  std::uint64_t internal_element_id;
  std::uint64_t gauss_point_index;
  std::uint64_t section_point_index;
  double x1;
  double x2;
  const char* source;
  double s11;
};

struct DiagnosticWriteRow {
  const char* severity;
  const char* code;
  const char* file;
  std::uint64_t line;
  const char* keyword;
  const char* entity_identity;
  const char* message;
};

bool IsOrthonormalRightHanded(
    const std::array<std::array<double, 3>, 3>& frame) {
  constexpr double kTolerance = 1.0e-12;
  const std::array<Vector3, 3> axes = {Vector3{frame[0U]}, Vector3{frame[1U]},
                                       Vector3{frame[2U]}};
  for (const Vector3& axis : axes) {
    if (!axis.IsFinite() || std::abs(axis.Dot(axis) - 1.0) > kTolerance) {
      return false;
    }
  }
  if (std::abs(axes[0U].Dot(axes[1U])) > kTolerance ||
      std::abs(axes[0U].Dot(axes[2U])) > kTolerance ||
      std::abs(axes[1U].Dot(axes[2U])) > kTolerance) {
    return false;
  }
  const Vector3 cross = axes[0U].Cross(axes[1U]);
  const double handedness = cross.Dot(axes[2U]);
  return handedness > 0.0 && std::abs(handedness - 1.0) <= kTolerance;
}

bool ComputeLocalAxes(const Domain& domain,
                      const EulerBeam3DDefinition& element, AxisSet& axes) {
  if (element.node_indices[0U] >= domain.Nodes().size() ||
      element.node_indices[1U] >= domain.Nodes().size() ||
      element.section_index >= domain.Sections().Size()) {
    return false;
  }
  const auto& first = domain.Nodes()[element.node_indices[0U]].coordinates;
  const auto& second = domain.Nodes()[element.node_indices[1U]].coordinates;
  const Vector3 first_position{first};
  const Vector3 second_position{second};
  const Vector3 guide{domain.Sections()[element.section_index].first_axis};
  const Vector3 delta = second_position - first_position;
  const double length = delta.Norm();
  if (!first_position.IsFinite() || !second_position.IsFinite() ||
      !guide.IsFinite() || !delta.IsFinite() || !std::isfinite(length) ||
      !(length > 0.0)) {
    return false;
  }
  const Vector3 x = delta / length;
  const double projection = guide.Dot(x);
  const Vector3 y_trial = guide - projection * x;
  const double y_norm = y_trial.Norm();
  if (!y_trial.IsFinite() || !std::isfinite(y_norm) || !(y_norm > 0.0)) {
    return false;
  }
  const Vector3 y = y_trial / y_norm;
  const Vector3 z = x.Cross(y);
  axes = {x[0U], x[1U], x[2U], y[0U], y[1U], y[2U], z[0U], z[1U], z[2U]};
  return IsFinite(axes);
}

Status ValidateShellWriterInput(const Domain& domain,
                                const AnalysisState& state) {
  if (!domain.BeamElements().Empty()) {
    return OutputFailure(
        "invalid-result-identity",
        "Schema v0 does not combine B33 and FESA-MITC4 element inventories.");
  }
  for (std::size_t index = 0U; index < domain.LinearElasticMaterials().Size();
       ++index) {
    const auto& material = domain.LinearElasticMaterials()[index];
    if (material.name.empty() || !IsValidUtf8(material.name) ||
        !std::isfinite(material.youngs_modulus) ||
        !std::isfinite(material.poisson_ratio) ||
        !(material.youngs_modulus > 0.0)) {
      return OutputFailure("invalid-result-identity",
                           "Every shell material requires a UTF-8 name and "
                           "finite constitutive data.");
    }
  }
  for (std::size_t index = 0U; index < domain.ShellSections().Size(); ++index) {
    const auto& section = domain.ShellSections()[index];
    if (section.name.empty() || !IsValidUtf8(section.name) ||
        section.material_index >= domain.LinearElasticMaterials().Size() ||
        !std::isfinite(section.thickness) || !(section.thickness > 0.0)) {
      return OutputFailure("invalid-result-identity",
                           "Every shell section requires stable material "
                           "identity and positive finite thickness.");
    }
  }
  for (std::size_t index = 0U; index < domain.ShellElements().Size(); ++index) {
    const auto& element = domain.ShellElements()[index];
    if ((element.source_type != ShellSourceElementType::kS4 &&
         element.source_type != ShellSourceElementType::kS4r) ||
        element.source_id.source_label <= 0 ||
        element.source_id.source_label_text.empty() ||
        !IsValidUtf8(element.source_id.instance_name) ||
        !IsValidUtf8(element.source_id.source_label_text) ||
        element.node_indices.size() != kShellNodeCount ||
        element.material_index >= domain.LinearElasticMaterials().Size() ||
        element.section_index >= domain.ShellSections().Size() ||
        domain.ShellSections()[element.section_index].material_index !=
            element.material_index) {
      return OutputFailure("invalid-result-identity",
                           "Every shell element requires stable source, "
                           "material, and section identity.");
    }
    auto sorted_nodes = element.node_indices;
    std::sort(sorted_nodes.begin(), sorted_nodes.end());
    if (sorted_nodes.back() >= domain.Nodes().size() ||
        std::adjacent_find(sorted_nodes.begin(), sorted_nodes.end()) !=
            sorted_nodes.end()) {
      return OutputFailure(
          "invalid-result-identity",
          "Every shell element requires four distinct valid node identities.");
    }
  }
  if (domain.ShellNodeInitialFrames().size() != domain.Nodes().size()) {
    return OutputFailure("invalid-result-identity",
                         "Shell output requires one initial director/frame per "
                         "source-ordered node.");
  }
  for (std::size_t node = 0U; node < domain.ShellNodeInitialFrames().size();
       ++node) {
    const auto& source = domain.ShellNodeInitialFrames()[node];
    const std::array<std::array<double, 3>, 3> frame{
        source.tangent_a, source.tangent_b, source.director};
    if (source.node_index != node || !IsOrthonormalRightHanded(frame)) {
      return OutputFailure("invalid-result-identity",
                           "Shell initial frames must be finite, orthonormal, "
                           "right-handed, and node ordered.");
    }
  }

  std::size_t expected_rows = 0U;
  if (!SizeProductFits(domain.ShellElements().Size(), kShellLocationCount,
                       expected_rows) ||
      state.ShellResults().size() != expected_rows) {
    return OutputFailure("invalid-result-rows",
                         "Shell output requires exactly GP1 through GP4 for "
                         "every shell element.");
  }
  const double gauss = 1.0 / std::sqrt(3.0);
  const std::array<ShellMidsurfaceLocation, kShellLocationCount> locations{
      ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2,
      ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4};
  const std::array<std::array<double, 2>, kShellLocationCount> coordinates{
      {{-gauss, -gauss}, {gauss, -gauss}, {gauss, gauss}, {-gauss, gauss}}};
  const std::array<ShellSectionPosition, kShellSectionPositionCount> positions{
      ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle,
      ShellSectionPosition::kTop};
  constexpr std::array<double, kShellSectionPositionCount> kZeta{-1.0, 0.0,
                                                                 1.0};
  for (std::size_t row_index = 0U; row_index < state.ShellResults().size();
       ++row_index) {
    const auto& row = state.ShellResults()[row_index];
    const std::size_t element = row_index / kShellLocationCount;
    const std::size_t location = row_index % kShellLocationCount;
    if (row.element != element || row.location != locations[location] ||
        row.natural_coordinates != coordinates[location] ||
        !IsOrthonormalRightHanded(row.local_frame) ||
        !IsFinite(row.generalized_strain) || !IsFinite(row.section_resultant)) {
      return OutputFailure("invalid-result-rows",
                           "Shell result rows must be finite and preserve "
                           "element/GP/frame identity.");
    }
    for (std::size_t position = 0U; position < kShellSectionPositionCount;
         ++position) {
      if (row.stress[position].position != positions[position] ||
          row.stress[position].zeta != kZeta[position] ||
          !Vector3{row.stress[position].components}.IsFinite()) {
        return OutputFailure("invalid-result-rows",
                             "Shell stress rows must preserve BOTTOM, "
                             "MIDDLE, TOP identity.");
      }
    }
  }
  if (!std::isfinite(state.PhysicalStrainEnergy()) ||
      !IsFinite(state.Equilibrium()) ||
      !Vector3{state.VerificationMetrics()}.IsFinite()) {
    return OutputFailure(
        "invalid-result-rows",
        "Shell energy, equilibrium, and verification metrics must be finite.");
  }
  return Status::Ok();
}

Status ValidateWriterInput(const std::filesystem::path& output_path,
                           const Domain& domain, const AnalysisState& state,
                           const std::vector<Diagnostic>& diagnostics,
                           WriterModelData& model_data) {
  if (output_path.empty() || output_path.filename().empty()) {
    return OutputFailure("invalid-output-path",
                         "The HDF5 output path must name a file.");
  }
  if (state.Identity().step_name != kStepName ||
      state.Identity().frame_index != kFrameIndex) {
    return OutputFailure(
        "invalid-result-state",
        "Schema v0 requires literal Step-1 and frame index 0.");
  }

  const bool shell = IsShellDomain(domain);
  if (shell) {
    const Status shell_validation = ValidateShellWriterInput(domain, state);
    if (!shell_validation.IsOk()) {
      return shell_validation;
    }
  } else if (!state.ShellResults().empty()) {
    return OutputFailure("invalid-result-rows",
                         "Beam output cannot contain shell recovery rows.");
  }

  std::size_t full_dof_count = 0U;
  if (!SizeProductFits(domain.Nodes().size(), kDofsPerNode, full_dof_count)) {
    return OutputFailure("invalid-result-state",
                         "The nodal result shape overflows size_t.");
  }
  const std::array<const Vector*, 5> vectors = {
      &state.Displacement(), &state.ExternalForce(), &state.InternalForce(),
      &state.Residual(), &state.Reaction()};
  for (const Vector* vector : vectors) {
    if (vector->Size() != full_dof_count) {
      return OutputFailure(
          "invalid-result-state",
          "Every V0 analysis vector must have node_count*6 values.");
    }
    for (std::size_t index = 0U; index < vector->Size(); ++index) {
      if (!std::isfinite((*vector)[index])) {
        return OutputFailure("invalid-result-state",
                             "Every V0 analysis vector value must be finite.");
      }
    }
  }

  if (domain.SourcePath().empty() || domain.SourceContentIdentity().empty() ||
      !IsValidUtf8(domain.SourceContentIdentity())) {
    return OutputFailure(
        "invalid-result-identity",
        "Source path and UTF-8 content identity are required.");
  }
  for (const Node& node : domain.Nodes()) {
    if (node.source_id.source_label <= 0 ||
        node.source_id.source_label_text.empty() ||
        !IsValidUtf8(node.source_id.instance_name) ||
        !IsValidUtf8(node.source_id.source_label_text) ||
        !Vector3{node.coordinates}.IsFinite()) {
      return OutputFailure(
          "invalid-result-identity",
          "Every node requires finite coordinates and UTF-8 source identity.");
    }
  }

  model_data.beam_local_axes.clear();
  model_data.beam_local_axes.reserve(domain.BeamElements().Size());
  for (std::size_t index = 0U; index < domain.BeamElements().Size(); ++index) {
    const EulerBeam3DDefinition& element = domain.BeamElements()[index];
    AxisSet axes{};
    if (element.source_id.source_label <= 0 ||
        element.source_id.source_label_text.empty() ||
        !IsValidUtf8(element.source_id.instance_name) ||
        !IsValidUtf8(element.source_id.source_label_text) ||
        element.node_indices[0U] == element.node_indices[1U] ||
        element.material_index >= domain.LinearElasticMaterials().Size() ||
        !ComputeLocalAxes(domain, element, axes)) {
      return OutputFailure("invalid-result-identity",
                           "Every element requires valid source, connectivity, "
                           "property, and local-axis identity.");
    }
    model_data.beam_local_axes.push_back(axes);
  }

  std::size_t endpoint_count = 0U;
  std::size_t gauss_count = 0U;
  if (!SizeProductFits(domain.BeamElements().Size(), kEndpointCount,
                       endpoint_count) ||
      !SizeProductFits(domain.BeamElements().Size(), kGaussPointCount,
                       gauss_count) ||
      state.EndpointResults().size() != endpoint_count ||
      state.GaussResults().size() != gauss_count) {
    return OutputFailure(
        "invalid-result-rows",
        "Endpoint and Gauss row counts must match every element and location.");
  }
  for (std::size_t row_index = 0U; row_index < state.EndpointResults().size();
       ++row_index) {
    const EntityIndex expected_element =
        static_cast<EntityIndex>(row_index / kEndpointCount);
    const int expected_endpoint = static_cast<int>(row_index % kEndpointCount);
    const EndpointResultRow& row = state.EndpointResults()[row_index];
    const auto& element = domain.BeamElements()[expected_element];
    const auto& expected_node =
        domain.Nodes()[element.node_indices[static_cast<std::size_t>(
            expected_endpoint)]];
    if (row.element != expected_element || row.endpoint != expected_endpoint ||
        !SameIdentity(row.node, expected_node.source_id) ||
        !IsFinite(row.end_action) || !IsFinite(row.section_resultant)) {
      return OutputFailure("invalid-result-rows",
                           "Endpoint result rows must follow "
                           "element/endpoint order and identity.");
    }
  }
  for (std::size_t row_index = 0U; row_index < state.GaussResults().size();
       ++row_index) {
    const EntityIndex expected_element =
        static_cast<EntityIndex>(row_index / kGaussPointCount);
    const int expected_gauss_point =
        static_cast<int>(row_index % kGaussPointCount) + 1;
    const GaussResultRow& row = state.GaussResults()[row_index];
    if (row.element != expected_element ||
        row.gauss_point != expected_gauss_point ||
        !IsFinite(row.generalized_strain) ||
        !IsFinite(row.generalized_resultant)) {
      return OutputFailure(
          "invalid-result-rows",
          "Gauss result rows must follow element/Gauss order and identity.");
    }
  }

  std::size_t stress_index = 0U;
  for (std::size_t element_index = 0U;
       element_index < domain.BeamElements().Size(); ++element_index) {
    const auto& element = domain.BeamElements()[element_index];
    const auto& section_points =
        domain.Sections()[element.section_index].section_points;
    for (std::size_t gauss = 0U; gauss < kGaussPointCount; ++gauss) {
      const std::size_t count =
          section_points.empty() ? 1U : section_points.size();
      for (std::size_t point = 0U; point < count; ++point) {
        if (stress_index >= state.StressResults().size()) {
          return OutputFailure("invalid-result-rows",
                               "Axial stress rows are missing required "
                               "element/Gauss/section locations.");
        }
        const StressS11Row& row = state.StressResults()[stress_index++];
        const std::size_t expected_point =
            section_points.empty() ? 0U : point + 1U;
        const double expected_x1 =
            section_points.empty() ? 0.0 : section_points[point][0U];
        const double expected_x2 =
            section_points.empty() ? 0.0 : section_points[point][1U];
        const char* expected_source =
            section_points.empty() ? "fesa-default" : "input";
        if (row.element != static_cast<EntityIndex>(element_index) ||
            row.gauss_point != static_cast<int>(gauss + 1U) ||
            row.section_point != expected_point || row.x1 != expected_x1 ||
            row.x2 != expected_x2 || row.source != expected_source ||
            !IsValidUtf8(row.source) || !std::isfinite(row.x1) ||
            !std::isfinite(row.x2) || !std::isfinite(row.s11)) {
          return OutputFailure("invalid-result-rows",
                               "Axial stress rows must follow exact "
                               "element/Gauss/section identity.");
        }
      }
    }
  }
  if (stress_index != state.StressResults().size()) {
    return OutputFailure("invalid-result-rows",
                         "Axial stress output contains extra rows.");
  }

  for (const Diagnostic& diagnostic : diagnostics) {
    if (!IsValidUtf8(diagnostic.code) || !IsValidUtf8(diagnostic.keyword) ||
        !IsValidUtf8(diagnostic.entity_identity) ||
        !IsValidUtf8(diagnostic.message)) {
      return OutputFailure("invalid-result-diagnostic",
                           "Diagnostic text must be valid UTF-8.");
    }
  }

  auto analysis_model_result = AnalysisModel::Create(domain);
  if (!analysis_model_result.HasValue()) {
    return OutputFailure(
        "invalid-result-state",
        "The HDF5 writer could not reconstruct the active model view.");
  }
  const AnalysisModel analysis_model = std::move(analysis_model_result.Value());
  auto dof_result = DofManager::Create(analysis_model);
  if (!dof_result.HasValue()) {
    return OutputFailure(
        "invalid-result-state",
        "The HDF5 writer could not reconstruct stable constraint identity.");
  }
  const DofManager dofs = std::move(dof_result.Value());
  model_data.constraint_mask.assign(full_dof_count, 0U);
  model_data.prescribed_displacement.assign(full_dof_count, 0.0);
  if (dofs.ConstrainedDofs().size() != dofs.PrescribedValues().Size()) {
    return OutputFailure(
        "invalid-result-state",
        "Constraint identities and prescribed values have inconsistent sizes.");
  }
  for (std::size_t index = 0U; index < dofs.ConstrainedDofs().size(); ++index) {
    const std::size_t full_dof = dofs.ConstrainedDofs()[index];
    const double prescribed = dofs.PrescribedValues()[index];
    if (full_dof >= full_dof_count || !std::isfinite(prescribed)) {
      return OutputFailure("invalid-result-state",
                           "Constraint identity and prescribed values must "
                           "be finite and in range.");
    }
    model_data.constraint_mask[full_dof] = 1U;
    model_data.prescribed_displacement[full_dof] = prescribed;
  }
  return Status::Ok();
}

std::vector<double> FlattenEndpointValues(
    const std::vector<EndpointResultRow>& rows, const bool section_resultants) {
  const std::size_t components = section_resultants ? kGeneralizedComponentCount
                                                    : kEndActionComponentCount;
  std::vector<double> values;
  values.reserve(rows.size() * components);
  for (const auto& row : rows) {
    if (section_resultants) {
      values.insert(values.end(), row.section_resultant.begin(),
                    row.section_resultant.end());
    } else {
      values.insert(values.end(), row.end_action.begin(), row.end_action.end());
    }
  }
  return values;
}

std::vector<double> FlattenGaussValues(const std::vector<GaussResultRow>& rows,
                                       const bool resultants) {
  std::vector<double> values;
  values.reserve(rows.size() * kGeneralizedComponentCount);
  for (const auto& row : rows) {
    const auto& row_values =
        resultants ? row.generalized_resultant : row.generalized_strain;
    values.insert(values.end(), row_values.begin(), row_values.end());
  }
  return values;
}

void WriteStress(const hid_t file, const AnalysisState& state) {
  std::vector<StressWriteRow> rows;
  rows.reserve(state.StressResults().size());
  for (const auto& row : state.StressResults()) {
    rows.push_back({static_cast<std::uint64_t>(row.element),
                    static_cast<std::uint64_t>(row.gauss_point),
                    static_cast<std::uint64_t>(row.section_point), row.x1,
                    row.x2, row.source.c_str(), row.s11});
  }
  auto string_type = MakeUtf8StringType();
  Hdf5Handle file_type{
      RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(StressWriteRow)),
                    "Unable to create the stress file type."),
      H5Tclose};
  Hdf5Handle memory_type{
      RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(StressWriteRow)),
                    "Unable to create the stress memory type."),
      H5Tclose};
  const auto insert_fields = [&](const hid_t type, const hid_t integer_type,
                                 const hid_t double_type) {
    RequireHdf5(
        H5Tinsert(type, "internal_element_id",
                  HOFFSET(StressWriteRow, internal_element_id), integer_type),
        "Unable to define the stress element field.");
    RequireHdf5(
        H5Tinsert(type, "gauss_point_index",
                  HOFFSET(StressWriteRow, gauss_point_index), integer_type),
        "Unable to define the stress Gauss field.");
    RequireHdf5(
        H5Tinsert(type, "section_point_index",
                  HOFFSET(StressWriteRow, section_point_index), integer_type),
        "Unable to define the stress section-point field.");
    RequireHdf5(H5Tinsert(type, "x1", HOFFSET(StressWriteRow, x1), double_type),
                "Unable to define the stress x1 field.");
    RequireHdf5(H5Tinsert(type, "x2", HOFFSET(StressWriteRow, x2), double_type),
                "Unable to define the stress x2 field.");
    RequireHdf5(H5Tinsert(type, "source", HOFFSET(StressWriteRow, source),
                          string_type.Get()),
                "Unable to define the stress source field.");
    RequireHdf5(
        H5Tinsert(type, "S11", HOFFSET(StressWriteRow, s11), double_type),
        "Unable to define the stress S11 field.");
  };
  insert_fields(file_type.Get(), H5T_STD_U64LE, H5T_IEEE_F64LE);
  insert_fields(memory_type.Get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE);
  auto dataset = WriteCompoundDataset(
      file, "/steps/Step-1/frames/0/element/stress_s11", rows.size(),
      file_type.Get(), memory_type.Get(), rows.data());
  WriteResultAttributes(dataset.Get(), "S11", "force/length^2", "beam-local",
                        "section-point");
}

void WriteShellResultIdentity(const hid_t file, const std::string& path,
                              const bool uses_section_positions) {
  Hdf5Handle dataset{RequireHdf5Id(H5Dopen2(file, path.c_str(), H5P_DEFAULT),
                                   "Unable to reopen a shell result dataset."),
                     H5Dclose};
  WriteStringAttribute(dataset.Get(), "source_element_type_dataset",
                       "/model/elements.source_element_type");
  WriteStringAttribute(dataset.Get(), "internal_formulation", "FESA-MITC4");
  WriteStringAttribute(dataset.Get(), "midsurface_location_dataset",
                       "/model/shell/midsurface_locations");
  WriteStringAttribute(dataset.Get(), "local_frame_dataset",
                       "/steps/Step-1/frames/0/element/shell/local_frame");
  if (uses_section_positions) {
    WriteStringAttribute(dataset.Get(), "section_position_dataset",
                         "/model/shell/section_positions");
  }
}

void WriteShellResultDatasets(const hid_t file, const Domain& domain,
                              const AnalysisState& state) {
  std::vector<double> local_frames;
  std::vector<double> generalized_strains;
  std::vector<double> section_resultants;
  std::vector<double> stresses;
  local_frames.reserve(state.ShellResults().size() * 9U);
  generalized_strains.reserve(state.ShellResults().size() *
                              kShellGeneralizedComponentCount);
  section_resultants.reserve(state.ShellResults().size() *
                             kShellGeneralizedComponentCount);
  stresses.reserve(state.ShellResults().size() * kShellSectionPositionCount *
                   kShellStressComponentCount);
  for (const auto& row : state.ShellResults()) {
    for (const auto& axis : row.local_frame) {
      local_frames.insert(local_frames.end(), axis.begin(), axis.end());
    }
    generalized_strains.insert(generalized_strains.end(),
                               row.generalized_strain.begin(),
                               row.generalized_strain.end());
    section_resultants.insert(section_resultants.end(),
                              row.section_resultant.begin(),
                              row.section_resultant.end());
    for (const auto& position : row.stress) {
      stresses.insert(stresses.end(), position.components.begin(),
                      position.components.end());
    }
  }

  const hsize_t element_count =
      static_cast<hsize_t>(domain.ShellElements().Size());
  const std::string root = std::string{kStepRoot} + "/element/shell";
  const std::string frame_path = root + "/local_frame";
  WriteDoubleDataset(file, frame_path, {element_count, 4U, 3U, 3U},
                     local_frames.data(), local_frames.size(), "X,Y,Z", "1,1,1",
                     "global-cartesian", "shell-local-frame");
  {
    Hdf5Handle dataset{
        RequireHdf5Id(H5Dopen2(file, frame_path.c_str(), H5P_DEFAULT),
                      "Unable to reopen shell local-frame results."),
        H5Dclose};
    WriteStringAttribute(dataset.Get(), "axis_names", "E1,E2,E3");
    WriteStringAttribute(dataset.Get(), "source_element_type_dataset",
                         "/model/elements.source_element_type");
    WriteStringAttribute(dataset.Get(), "internal_formulation", "FESA-MITC4");
    WriteStringAttribute(dataset.Get(), "midsurface_location_dataset",
                         "/model/shell/midsurface_locations");
  }

  const std::string strain_path = root + "/generalized_strain";
  WriteDoubleDataset(
      file, strain_path, {element_count, 4U, 8U}, generalized_strains.data(),
      generalized_strains.size(), "E11,E22,G12,K11,K22,K12,G13,G23",
      "1,1,1,1/length,1/length,1/length,1,1", "shell-local", "midsurface");
  WriteShellResultIdentity(file, strain_path, false);

  const std::string resultant_path = root + "/section_resultant";
  WriteDoubleDataset(file, resultant_path, {element_count, 4U, 8U},
                     section_resultants.data(), section_resultants.size(),
                     "N11,N22,N12,M11,M22,M12,Q13,Q23",
                     "force/length,force/length,force/"
                     "length,force,force,force,force/length,force/length",
                     "shell-local", "midsurface");
  WriteShellResultIdentity(file, resultant_path, false);

  const std::string stress_path = root + "/stress";
  WriteDoubleDataset(file, stress_path, {element_count, 4U, 3U, 3U},
                     stresses.data(), stresses.size(), "S11,S22,S12",
                     "force/length^2,force/length^2,force/length^2",
                     "shell-local", "section-position");
  WriteShellResultIdentity(file, stress_path, true);

  const double energy = state.PhysicalStrainEnergy();
  WriteDoubleDataset(file, std::string{kStepRoot} + "/global/energy", {1U},
                     &energy, 1U, "PHYSICAL_STRAIN_ENERGY", "force*length",
                     "global", "global");
  WriteDoubleDataset(file, std::string{kStepRoot} + "/global/equilibrium", {6U},
                     state.Equilibrium().data(), state.Equilibrium().size(),
                     "FORCE_1,FORCE_2,FORCE_3,MOMENT_1,MOMENT_2,MOMENT_3",
                     "force,force,force,force*length,force*length,force*length",
                     "global-cartesian", "global-origin");
  const std::string metrics_path =
      std::string{kStepRoot} + "/global/verification_metrics";
  WriteDoubleDataset(file, metrics_path, {3U},
                     state.VerificationMetrics().data(),
                     state.VerificationMetrics().size(),
                     "FREE_RESIDUAL_NORMALIZED,FORCE_BALANCE_NORMALIZED,"
                     "MOMENT_BALANCE_NORMALIZED",
                     "1,1,1", "global", "verification");
  {
    Hdf5Handle dataset{
        RequireHdf5Id(H5Dopen2(file, metrics_path.c_str(), H5P_DEFAULT),
                      "Unable to reopen shell verification metrics."),
        H5Dclose};
    WriteStringAttribute(
        dataset.Get(), "metric_definition_ids",
        "free-residual-l2-over-max-free-force-l2,force-balance-l2-over-max-"
        "force-sum,moment-balance-l2-over-max-moment-sum");
    WriteStringAttribute(dataset.Get(), "acceptance_thresholds",
                         "1e-10,1e-10,1e-10");
  }
}

void WriteDiagnostics(const hid_t file,
                      const std::vector<Diagnostic>& input_diagnostics) {
  std::vector<Diagnostic> diagnostics = input_diagnostics;
  SortDiagnostics(diagnostics);
  std::vector<std::string> files;
  files.reserve(diagnostics.size());
  for (const auto& diagnostic : diagnostics) {
    files.push_back(NormalizedPathString(diagnostic.location.file));
  }
  std::vector<DiagnosticWriteRow> rows;
  rows.reserve(diagnostics.size());
  for (std::size_t index = 0U; index < diagnostics.size(); ++index) {
    const auto& diagnostic = diagnostics[index];
    rows.push_back(
        {diagnostic.severity == Severity::kWarning ? "warning" : "error",
         diagnostic.code.c_str(), files[index].c_str(),
         static_cast<std::uint64_t>(diagnostic.location.line),
         diagnostic.keyword.c_str(), diagnostic.entity_identity.c_str(),
         diagnostic.message.c_str()});
  }

  auto string_type = MakeUtf8StringType();
  Hdf5Handle file_type{
      RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(DiagnosticWriteRow)),
                    "Unable to create the diagnostic file type."),
      H5Tclose};
  Hdf5Handle memory_type{
      RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(DiagnosticWriteRow)),
                    "Unable to create the diagnostic memory type."),
      H5Tclose};
  const auto insert_fields = [&](const hid_t type, const hid_t integer_type) {
    RequireHdf5(
        H5Tinsert(type, "severity", HOFFSET(DiagnosticWriteRow, severity),
                  string_type.Get()),
        "Unable to define diagnostic severity.");
    RequireHdf5(H5Tinsert(type, "code", HOFFSET(DiagnosticWriteRow, code),
                          string_type.Get()),
                "Unable to define diagnostic code.");
    RequireHdf5(H5Tinsert(type, "file", HOFFSET(DiagnosticWriteRow, file),
                          string_type.Get()),
                "Unable to define diagnostic file.");
    RequireHdf5(H5Tinsert(type, "line", HOFFSET(DiagnosticWriteRow, line),
                          integer_type),
                "Unable to define diagnostic line.");
    RequireHdf5(H5Tinsert(type, "keyword", HOFFSET(DiagnosticWriteRow, keyword),
                          string_type.Get()),
                "Unable to define diagnostic keyword.");
    RequireHdf5(H5Tinsert(type, "entity_identity",
                          HOFFSET(DiagnosticWriteRow, entity_identity),
                          string_type.Get()),
                "Unable to define diagnostic entity identity.");
    RequireHdf5(H5Tinsert(type, "message", HOFFSET(DiagnosticWriteRow, message),
                          string_type.Get()),
                "Unable to define diagnostic message.");
  };
  insert_fields(file_type.Get(), H5T_STD_U64LE);
  insert_fields(memory_type.Get(), H5T_NATIVE_UINT64);
  (void)WriteCompoundDataset(file, "/diagnostics", rows.size(), file_type.Get(),
                             memory_type.Get(), rows.data());
}

void WriteResultDatasets(const hid_t file, const Domain& domain,
                         const AnalysisState& state) {
  const std::vector<hsize_t> nodal_dimensions = {
      static_cast<hsize_t>(domain.Nodes().size()), kDofsPerNode};
  WriteDoubleDataset(file, std::string{kStepRoot} + "/nodal/displacement",
                     nodal_dimensions, state.Displacement().Data(),
                     state.Displacement().Size(), "UX,UY,UZ,URX,URY,URZ",
                     "length,length,length,radian,radian,radian",
                     "global-cartesian", "nodal");
  WriteDoubleDataset(file, std::string{kStepRoot} + "/nodal/reaction",
                     nodal_dimensions, state.Reaction().Data(),
                     state.Reaction().Size(), "RF1,RF2,RF3,RM1,RM2,RM3",
                     "force,force,force,force*length,force*length,force*length",
                     "global-cartesian", "nodal");

  if (IsShellDomain(domain)) {
    WriteShellResultDatasets(file, domain, state);
    return;
  }

  const std::vector<hsize_t> endpoint_action_dimensions = {
      static_cast<hsize_t>(domain.BeamElements().Size()), kEndpointCount,
      kEndActionComponentCount};
  const std::vector<hsize_t> generalized_dimensions = {
      static_cast<hsize_t>(domain.BeamElements().Size()), kEndpointCount,
      kGeneralizedComponentCount};
  const auto end_actions =
      FlattenEndpointValues(state.EndpointResults(), false);
  WriteDoubleDataset(file, std::string{kStepRoot} + "/element/end_force_local",
                     endpoint_action_dimensions, end_actions.data(),
                     end_actions.size(), "FX,FY,FZ,MX,MY,MZ",
                     "force,force,force,force*length,force*length,force*length",
                     "beam-local", "endpoint-outward-action");
  const auto section_resultants =
      FlattenEndpointValues(state.EndpointResults(), true);
  WriteDoubleDataset(file,
                     std::string{kStepRoot} + "/element/section_resultant",
                     generalized_dimensions, section_resultants.data(),
                     section_resultants.size(), "N,T,My,Mz",
                     "force,force*length,force*length,force*length",
                     "beam-local", "endpoint-positive-local-x-section-cut");
  const auto generalized_strains =
      FlattenGaussValues(state.GaussResults(), false);
  WriteDoubleDataset(
      file, std::string{kStepRoot} + "/element/generalized_strain",
      generalized_dimensions, generalized_strains.data(),
      generalized_strains.size(), "epsilon0,kappa_x,kappa_y,kappa_z",
      "1,1/length,1/length,1/length", "beam-local", "integration-point");
  const auto generalized_resultants =
      FlattenGaussValues(state.GaussResults(), true);
  WriteDoubleDataset(file,
                     std::string{kStepRoot} + "/element/generalized_resultant",
                     generalized_dimensions, generalized_resultants.data(),
                     generalized_resultants.size(), "N,T,My,Mz",
                     "force,force*length,force*length,force*length",
                     "beam-local", "integration-point");
  WriteStress(file, state);
}

std::vector<double> CopyVector(const Vector& source) {
  std::vector<double> values(source.Size());
  for (std::size_t index = 0U; index < source.Size(); ++index) {
    values[index] = source[index];
  }
  return values;
}

std::vector<DoubleDatasetPlan> BuildResultDatasetPlans(
    const Domain& domain, const AnalysisState& state) {
  std::vector<DoubleDatasetPlan> datasets;
  const std::vector<std::size_t> nodal_dimensions{domain.Nodes().size(),
                                                  kDofsPerNode};
  datasets.push_back({std::string{kStepRoot} + "/nodal/displacement",
                      nodal_dimensions, CopyVector(state.Displacement()),
                      "UX,UY,UZ,URX,URY,URZ",
                      "length,length,length,radian,radian,radian",
                      "global-cartesian", "nodal"});
  datasets.push_back(
      {std::string{kStepRoot} + "/nodal/reaction", nodal_dimensions,
       CopyVector(state.Reaction()), "RF1,RF2,RF3,RM1,RM2,RM3",
       "force,force,force,force*length,force*length,force*length",
       "global-cartesian", "nodal"});
  if (IsShellDomain(domain)) {
    std::vector<double> local_frames;
    std::vector<double> generalized_strains;
    std::vector<double> section_resultants;
    std::vector<double> stresses;
    for (const auto& row : state.ShellResults()) {
      for (const auto& axis : row.local_frame) {
        local_frames.insert(local_frames.end(), axis.begin(), axis.end());
      }
      generalized_strains.insert(generalized_strains.end(),
                                 row.generalized_strain.begin(),
                                 row.generalized_strain.end());
      section_resultants.insert(section_resultants.end(),
                                row.section_resultant.begin(),
                                row.section_resultant.end());
      for (const auto& position : row.stress) {
        stresses.insert(stresses.end(), position.components.begin(),
                        position.components.end());
      }
    }
    const std::size_t element_count = domain.ShellElements().Size();
    const std::string root = std::string{kStepRoot} + "/element/shell";
    datasets.push_back({root + "/local_frame",
                        {element_count, 4U, 3U, 3U},
                        std::move(local_frames),
                        "X,Y,Z",
                        "1,1,1",
                        "global-cartesian",
                        "shell-local-frame"});
    datasets.push_back({root + "/generalized_strain",
                        {element_count, 4U, 8U},
                        std::move(generalized_strains),
                        "E11,E22,G12,K11,K22,K12,G13,G23",
                        "1,1,1,1/length,1/length,1/length,1,1",
                        "shell-local",
                        "midsurface"});
    datasets.push_back({root + "/section_resultant",
                        {element_count, 4U, 8U},
                        std::move(section_resultants),
                        "N11,N22,N12,M11,M22,M12,Q13,Q23",
                        "force/length,force/length,force/length,force,force,"
                        "force,force/length,force/length",
                        "shell-local",
                        "midsurface"});
    datasets.push_back({root + "/stress",
                        {element_count, 4U, 3U, 3U},
                        std::move(stresses),
                        "S11,S22,S12",
                        "force/length^2,force/length^2,force/length^2",
                        "shell-local",
                        "section-position"});
    datasets.push_back({std::string{kStepRoot} + "/global/energy",
                        {1U},
                        {state.PhysicalStrainEnergy()},
                        "PHYSICAL_STRAIN_ENERGY",
                        "force*length",
                        "global",
                        "global"});
    datasets.push_back({std::string{kStepRoot} + "/global/equilibrium",
                        {6U},
                        std::vector<double>(state.Equilibrium().begin(),
                                            state.Equilibrium().end()),
                        "FORCE_1,FORCE_2,FORCE_3,MOMENT_1,MOMENT_2,MOMENT_3",
                        "force,force,force,force*length,force*length,"
                        "force*length",
                        "global-cartesian",
                        "global-origin"});
    datasets.push_back(
        {std::string{kStepRoot} + "/global/verification_metrics",
         {3U},
         std::vector<double>(state.VerificationMetrics().begin(),
                             state.VerificationMetrics().end()),
         "FREE_RESIDUAL_NORMALIZED,FORCE_BALANCE_NORMALIZED,MOMENT_BALANCE_"
         "NORMALIZED",
         "1,1,1",
         "global",
         "verification"});
    return datasets;
  }

  const std::size_t element_count = domain.BeamElements().Size();
  datasets.push_back({std::string{kStepRoot} + "/element/end_force_local",
                      {element_count, kEndpointCount, kEndActionComponentCount},
                      FlattenEndpointValues(state.EndpointResults(), false),
                      "FX,FY,FZ,MX,MY,MZ",
                      "force,force,force,force*length,force*length,"
                      "force*length",
                      "beam-local",
                      "endpoint-outward-action"});
  datasets.push_back(
      {std::string{kStepRoot} + "/element/section_resultant",
       {element_count, kEndpointCount, kGeneralizedComponentCount},
       FlattenEndpointValues(state.EndpointResults(), true),
       "N,T,My,Mz",
       "force,force*length,force*length,force*length",
       "beam-local",
       "endpoint-positive-local-x-section-cut"});
  datasets.push_back(
      {std::string{kStepRoot} + "/element/generalized_strain",
       {element_count, kGaussPointCount, kGeneralizedComponentCount},
       FlattenGaussValues(state.GaussResults(), false),
       "epsilon0,kappa_x,kappa_y,kappa_z",
       "1,1/length,1/length,1/length",
       "beam-local",
       "integration-point"});
  datasets.push_back(
      {std::string{kStepRoot} + "/element/generalized_resultant",
       {element_count, kGaussPointCount, kGeneralizedComponentCount},
       FlattenGaussValues(state.GaussResults(), true),
       "N,T,My,Mz",
       "force,force*length,force*length,force*length",
       "beam-local",
       "integration-point"});
  return datasets;
}

Status ValidateFiniteInventory(const std::vector<DoubleDatasetPlan>& datasets) {
  std::set<std::string> paths;
  for (const DoubleDatasetPlan& dataset : datasets) {
    std::size_t expected_values = 1U;
    bool valid_shape = !dataset.path.empty() && !dataset.dimensions.empty();
    for (const std::size_t dimension : dataset.dimensions) {
      std::size_t next = 0U;
      if (!SizeProductFits(expected_values, dimension, next)) {
        valid_shape = false;
        break;
      }
      expected_values = next;
    }
    const bool finite =
        std::all_of(dataset.values.begin(), dataset.values.end(),
                    [](const double value) { return std::isfinite(value); });
    if (!valid_shape || expected_values != dataset.values.size() || !finite ||
        !paths.insert(dataset.path).second) {
      return OutputFailure(
          "invalid-result-dataset-plan",
          "Result dataset plans require unique paths, exact shapes, and "
          "finite values.");
    }
  }
  return Status::Ok();
}

Result<WriterPlan> BuildWriterPlan(const std::filesystem::path& output_path,
                                   const Domain& domain,
                                   const AnalysisState& state,
                                   const std::vector<Diagnostic>& diagnostics) {
  WriterPlan plan;
  const Status validation = ValidateWriterInput(output_path, domain, state,
                                                diagnostics, plan.model_data);
  if (!validation.IsOk()) {
    return Result<WriterPlan>::Failure(validation);
  }
  plan.result_datasets = BuildResultDatasetPlans(domain, state);
  const Status inventory = ValidateFiniteInventory(plan.result_datasets);
  if (!inventory.IsOk()) {
    return Result<WriterPlan>::Failure(inventory);
  }
  return Result<WriterPlan>::Success(std::move(plan));
}

void WriteResults(const hid_t file, const Domain& domain,
                  const AnalysisState& state,
                  const std::vector<Diagnostic>& diagnostics,
                  const WriterPlan& plan) {
  const Status inventory = ValidateFiniteInventory(plan.result_datasets);
  if (!inventory.IsOk()) {
    throw Hdf5Failure{inventory.Diagnostics().front().message};
  }
  WriteResultDatasets(file, domain, state);
  WriteDiagnostics(file, diagnostics);
}

}  // namespace fesa::hdf5_internal