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

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

#include "fesa/analysis/analysis_model.h"
#include "fesa/build_info.h"
#include "fesa/elements/finite_rotation.h"
#include "fesa/fem/dof_manager.h"
#include "fesa/math/vector3.h"
#include "io/hdf5/hdf5_atomic_file.h"
#include "io/hdf5/hdf5_primitives.h"
#include "io/hdf5/hdf5_raii.h"
#include "io/hdf5/hdf5_self_check.h"

namespace fesa {
namespace {

using hdf5_internal::Hdf5Failure;
using hdf5_internal::Hdf5Handle;

constexpr std::size_t kDofsPerNode = 6U;
constexpr std::size_t kShellLocationsPerElement = 4U;
constexpr std::size_t kSectionPositions = 3U;

struct NonlinearModelData {
  std::vector<std::uint8_t> constraint_mask;
  std::vector<double> prescribed_start;
  std::vector<double> prescribed_target;
};

struct NodeWriteRow {
  std::uint64_t internal_node_id;
  const char* instance_name;
  const char* source_label;
  double coordinates[3];
};

struct ElementWriteRow {
  std::uint64_t internal_element_id;
  const char* instance_name;
  const char* source_label;
  const char* source_element_type;
  const char* internal_formulation;
  std::uint64_t node_internal_ids[4];
  std::uint64_t shell_section_internal_id;
  std::uint64_t material_internal_id;
};

struct MaterialWriteRow {
  std::uint64_t internal_material_id;
  const char* name;
  double youngs_modulus;
  double poisson_ratio;
};

struct SectionWriteRow {
  std::uint64_t internal_section_id;
  const char* source_file;
  std::uint64_t source_line;
  const char* source_elset;
  std::uint64_t material_internal_id;
  double thickness;
};

struct AcceptedFrameWriteRow {
  std::uint64_t frame_index;
  std::uint64_t increment_index;
  double load_factor;
  std::uint64_t iteration_count;
  std::uint8_t is_factor_one;
};

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;
};

template <std::size_t Size>
bool IsFinite(const std::array<double, Size>& values) {
  return std::all_of(values.begin(), values.end(),
                     [](const double value) { return std::isfinite(value); });
}

bool IsFinite(const Vector& values) {
  for (std::size_t index = 0U; index < values.Size(); ++index) {
    if (!std::isfinite(values[index])) {
      return false;
    }
  }
  return true;
}

bool SameIdentity(const SourceEntityId& left, const SourceEntityId& right) {
  return left.instance_name == right.instance_name &&
         left.source_label == right.source_label &&
         left.source_label_text == right.source_label_text;
}

bool IsFinite(const ShellResultRow& row) {
  if (!IsFinite(row.natural_coordinates) || !IsFinite(row.generalized_strain) ||
      !IsFinite(row.section_resultant)) {
    return false;
  }
  for (const auto& axis : row.local_frame) {
    if (!IsFinite(axis)) {
      return false;
    }
  }
  for (const auto& stress : row.stress) {
    if (!std::isfinite(stress.zeta) || !IsFinite(stress.components)) {
      return false;
    }
  }
  return true;
}

Status Failure(const char* code, const std::string& message) {
  return hdf5_internal::OutputFailure(code, message);
}

Status ValidateDomain(const Domain& domain) {
  if (!domain.BeamElements().Empty() || domain.ShellElements().Empty() ||
      domain.Steps().Size() != 1U || !domain.Steps()[0U].IsNonlinearStatic() ||
      domain.Nodes().empty() || domain.SourcePath().empty() ||
      domain.SourceContentIdentity().empty() ||
      !hdf5_internal::IsValidUtf8(domain.SourceContentIdentity())) {
    return Failure("invalid-result-identity",
                   "Schema v1 requires one nonlinear shell Domain identity.");
  }
  if (domain.ShellNodeInitialFrames().size() != domain.Nodes().size()) {
    return Failure("invalid-result-identity",
                   "Every nonlinear output node requires an initial frame.");
  }
  for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) {
    const Node& source = domain.Nodes()[node];
    const auto& frame = domain.ShellNodeInitialFrames()[node];
    if (source.source_id.source_label <= 0 ||
        source.source_id.source_label_text.empty() ||
        !hdf5_internal::IsValidUtf8(source.source_id.instance_name) ||
        !hdf5_internal::IsValidUtf8(source.source_id.source_label_text) ||
        !Vector3{source.coordinates}.IsFinite() ||
        frame.node_index != static_cast<EntityIndex>(node) ||
        !IsFinite(frame.tangent_a) || !IsFinite(frame.tangent_b) ||
        !IsFinite(frame.director)) {
      return Failure("invalid-result-identity",
                     "Nonlinear node identity or initial frame is invalid.");
    }
  }
  for (std::size_t element = 0U; element < domain.ShellElements().Size();
       ++element) {
    const auto& source = domain.ShellElements()[element];
    if (source.source_id.source_label <= 0 ||
        source.source_id.source_label_text.empty() ||
        !hdf5_internal::IsValidUtf8(source.source_id.instance_name) ||
        !hdf5_internal::IsValidUtf8(source.source_id.source_label_text) ||
        source.node_indices.size() != 4U ||
        source.section_index >= domain.ShellSections().Size() ||
        source.material_index >= domain.LinearElasticMaterials().Size() ||
        domain.ShellSections()[source.section_index].material_index !=
            source.material_index) {
      return Failure("invalid-result-identity",
                     "Nonlinear shell element identity is invalid.");
    }
    for (const EntityIndex node : source.node_indices) {
      if (node >= domain.Nodes().size()) {
        return Failure("invalid-result-identity",
                       "Nonlinear shell connectivity is out of range.");
      }
    }
  }
  return Status::Ok();
}

Status BuildModelData(const Domain& domain, NonlinearModelData& model_data) {
  auto model_result = AnalysisModel::Create(domain);
  if (!model_result.HasValue()) {
    return Failure("invalid-result-state",
                   "The nonlinear writer could not create an AnalysisModel.");
  }
  auto dofs_result = DofManager::Create(model_result.Value());
  if (!dofs_result.HasValue()) {
    return Failure("invalid-result-state",
                   "The nonlinear writer could not create a DofManager.");
  }
  const DofManager dofs = std::move(dofs_result.Value());
  std::size_t full_dof_count = 0U;
  if (!hdf5_internal::SizeProductFits(domain.Nodes().size(), kDofsPerNode,
                                      full_dof_count) ||
      dofs.FullDofCount() != full_dof_count ||
      dofs.ConstrainedDofs().size() != dofs.PrescribedValues().Size()) {
    return Failure("invalid-result-state",
                   "Nonlinear constraint identities have an invalid size.");
  }
  model_data.constraint_mask.assign(full_dof_count, 0U);
  model_data.prescribed_start.assign(full_dof_count, 0.0);
  model_data.prescribed_target.assign(full_dof_count, 0.0);
  for (std::size_t index = 0U; index < dofs.ConstrainedDofs().size(); ++index) {
    const std::size_t full_dof = dofs.ConstrainedDofs()[index];
    const double target = dofs.PrescribedValues()[index];
    if (full_dof >= full_dof_count || !std::isfinite(target)) {
      return Failure("invalid-result-state",
                     "Nonlinear prescribed values must be finite and in range.");
    }
    model_data.constraint_mask[full_dof] = 1U;
    model_data.prescribed_target[full_dof] = target;
  }
  return Status::Ok();
}

Status ValidateState(const Domain& domain, const NonlinearAnalysisState& state) {
  const auto& frames = state.AcceptedFrames();
  if (frames.empty()) {
    return Failure("invalid-result-state",
                   "Schema v1 requires at least one accepted frame.");
  }
  std::size_t full_dof_count = 0U;
  std::size_t shell_row_count = 0U;
  if (!hdf5_internal::SizeProductFits(domain.Nodes().size(), kDofsPerNode,
                                      full_dof_count) ||
      !hdf5_internal::SizeProductFits(domain.ShellElements().Size(),
                                      kShellLocationsPerElement,
                                      shell_row_count)) {
    return Failure("invalid-result-state",
                   "Nonlinear result dimensions overflow size_t.");
  }
  const double gauss = 1.0 / std::sqrt(3.0);
  const std::array<ShellMidsurfaceLocation, kShellLocationsPerElement>
      locations{ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2,
                ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4};
  const std::array<std::array<double, 2>, kShellLocationsPerElement>
      coordinates{{{-gauss, -gauss}, {gauss, -gauss},
                   {gauss, gauss}, {-gauss, gauss}}};
  for (std::size_t frame_index = 0U; frame_index < frames.size();
       ++frame_index) {
    const AcceptedFrameState& frame = frames[frame_index];
    const bool final_frame = frame_index + 1U == frames.size();
    if (frame.frame_index != frame_index || frame.increment_index == 0U ||
        !std::isfinite(frame.load_factor) || frame.load_factor <= 0.0 ||
        frame.load_factor > 1.0 || frame.iteration_count == 0U ||
        frame.is_factor_one != final_frame ||
        (final_frame && frame.load_factor != 1.0) ||
        (!final_frame && frame.load_factor >= 1.0) ||
        (frame_index > 0U &&
         (frame.increment_index <= frames[frame_index - 1U].increment_index ||
          frame.load_factor <= frames[frame_index - 1U].load_factor))) {
      return Failure("invalid-result-state",
                     "Accepted-frame ordering or factor-one identity is invalid.");
    }
    if (frame.nodal_internal_order.size() != domain.Nodes().size() ||
        frame.nodal_source_ids.size() != domain.Nodes().size() ||
        frame.nodal_triad.size() != domain.Nodes().size() ||
        frame.nodal_director.size() != domain.Nodes().size() ||
        frame.shell_results.size() != shell_row_count ||
        !std::isfinite(frame.physical_strain_energy) ||
        !IsFinite(frame.global_equilibrium)) {
      return Failure("invalid-result-state",
                     "Accepted-frame identity or physical inventory is invalid.");
    }
    const std::array<const Vector*, 7U> vectors{
        &frame.generalized_coordinates, &frame.displacement,
        &frame.rotation_vector,         &frame.internal_force,
        &frame.external_force,          &frame.residual,
        &frame.reaction};
    for (const Vector* values : vectors) {
      if (values->Size() != full_dof_count || !IsFinite(*values)) {
        return Failure("invalid-result-state",
                       "Accepted-frame vectors must be finite node-by-six data.");
      }
    }
    for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) {
      if (frame.nodal_internal_order[node] != node ||
          !SameIdentity(frame.nodal_source_ids[node],
                        domain.Nodes()[node].source_id) ||
          !IsFinite(frame.nodal_triad[node][0U]) ||
          !IsFinite(frame.nodal_triad[node][1U]) ||
          !IsFinite(frame.nodal_triad[node][2U]) ||
          !IsFinite(frame.nodal_director[node])) {
        return Failure("invalid-result-state",
                       "Accepted-frame nodal identity is invalid.");
      }
    }
    for (std::size_t row_index = 0U; row_index < shell_row_count;
         ++row_index) {
      const ShellResultRow& row = frame.shell_results[row_index];
      const std::size_t element = row_index / kShellLocationsPerElement;
      const std::size_t location = row_index % kShellLocationsPerElement;
      if (row.element != element || row.location != locations[location] ||
          row.natural_coordinates != coordinates[location] || !IsFinite(row)) {
        return Failure("invalid-result-state",
                       "Accepted shell rows must preserve element and GP order.");
      }
    }
  }
  return Status::Ok();
}

void WriteCommonNumericAttributes(const hid_t dataset,
                                  const std::string& component_names,
                                  const std::string& component_units,
                                  const std::string& coordinate_system,
                                  const std::string& location) {
  hdf5_internal::WriteUint64Attribute(dataset, "schema_version", 1U);
  hdf5_internal::WriteStringAttribute(dataset, "component_names",
                                      component_names);
  hdf5_internal::WriteStringAttribute(dataset, "component_unit_dimensions",
                                      component_units);
  hdf5_internal::WriteStringAttribute(dataset, "coordinate_system",
                                      coordinate_system);
  hdf5_internal::WriteStringAttribute(dataset, "location", location);
}

void WriteFrameAttributes(const hid_t dataset, const AcceptedFrameState& frame) {
  hdf5_internal::WriteStringAttribute(dataset, "step_name", "Step-1");
  hdf5_internal::WriteUint64Attribute(dataset, "frame_index", frame.frame_index);
  hdf5_internal::WriteUint64Attribute(dataset, "increment_index",
                                      frame.increment_index);
  hdf5_internal::WriteDoubleAttribute(dataset, "load_factor",
                                      frame.load_factor);
  hdf5_internal::WriteUint8Attribute(dataset, "is_factor_one",
                                     frame.is_factor_one ? 1U : 0U);
}

Hdf5Handle WriteModelDoubleDataset(
    const hid_t file, const std::string& path,
    const std::vector<hsize_t>& dimensions, const double* values,
    const std::size_t value_count, const std::string& component_names,
    const std::string& component_units, const std::string& coordinate_system,
    const std::string& location) {
  auto dataset = hdf5_internal::WriteDoubleValues(file, path, dimensions, values,
                                                  value_count);
  WriteCommonNumericAttributes(dataset.Get(), component_names, component_units,
                               coordinate_system, location);
  return dataset;
}

Hdf5Handle WriteFrameDoubleDataset(
    const hid_t file, const std::string& path,
    const std::vector<hsize_t>& dimensions, const double* values,
    const std::size_t value_count, const std::string& component_names,
    const std::string& component_units, const std::string& coordinate_system,
    const std::string& location, const AcceptedFrameState& frame) {
  auto dataset = WriteModelDoubleDataset(file, path, dimensions, values,
                                         value_count, component_names,
                                         component_units, coordinate_system,
                                         location);
  WriteFrameAttributes(dataset.Get(), frame);
  return dataset;
}

void WriteMetadata(const hid_t file, const Domain& domain) {
  auto metadata = hdf5_internal::CreateGroup(file, "/metadata");
  hdf5_internal::WriteUint64Attribute(metadata.Get(), "schema_version", 1U);
  hdf5_internal::WriteStringAttribute(metadata.Get(), "feature_id",
                                      "nonlinear-static-newton-raphson");
  hdf5_internal::WriteStringAttribute(metadata.Get(), "solver_version",
                                      std::string{SolverVersion()});
  hdf5_internal::WriteStringAttribute(metadata.Get(), "source_input_identity",
                                      hdf5_internal::SourceInputIdentity(domain));
  hdf5_internal::WriteStringAttribute(metadata.Get(), "unit_system_label",
                                      "user-consistent-unspecified");
  hdf5_internal::WriteStringAttribute(
      metadata.Get(), "coordinate_convention",
      "global-cartesian; shell-local=fixed-initial-reference; "
      "positive-thickness=+zeta");
  hdf5_internal::WriteStringAttribute(metadata.Get(), "analysis_procedure",
                                      "nonlinear-static-full-newton");
  hdf5_internal::WriteStringAttribute(metadata.Get(), "internal_formulation",
                                      "FESA-MITC4-total-lagrangian");
  hdf5_internal::WriteStringAttribute(metadata.Get(), "integration_rule",
                                      "2x2x2-gauss; mitc4-edge-midpoint-shear");
  hdf5_internal::WriteStringAttribute(metadata.Get(), "step_name", "Step-1");
  hdf5_internal::WriteStringAttribute(metadata.Get(), "frame_indexing",
                                      "accepted-increment-order-zero-based");
  hdf5_internal::WriteUint8Attribute(metadata.Get(), "factor_zero_is_frame",
                                     0U);
  hdf5_internal::WriteDoubleAttribute(metadata.Get(), "factor_zero_load_factor",
                                      0.0);
}

void WriteNodes(const hid_t file, const Domain& domain) {
  std::vector<NodeWriteRow> rows;
  rows.reserve(domain.Nodes().size());
  for (std::size_t index = 0U; index < domain.Nodes().size(); ++index) {
    const Node& node = domain.Nodes()[index];
    rows.push_back({static_cast<std::uint64_t>(index),
                    node.source_id.instance_name.c_str(),
                    node.source_id.source_label_text.c_str(),
                    {node.coordinates[0U], node.coordinates[1U],
                     node.coordinates[2U]}});
  }
  auto string_type = hdf5_internal::MakeUtf8StringType();
  const hsize_t coordinate_dimensions[] = {3U};
  Hdf5Handle file_coordinates{hdf5_internal::RequireHdf5Id(
                                  H5Tarray_create2(H5T_IEEE_F64LE, 1,
                                                   coordinate_dimensions),
                                  "Unable to create node coordinate type."),
                              H5Tclose};
  Hdf5Handle memory_coordinates{hdf5_internal::RequireHdf5Id(
                                    H5Tarray_create2(H5T_NATIVE_DOUBLE, 1,
                                                     coordinate_dimensions),
                                    "Unable to create node memory type."),
                                H5Tclose};
  Hdf5Handle file_type{hdf5_internal::RequireHdf5Id(
                           H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)),
                           "Unable to create nonlinear node file type."),
                       H5Tclose};
  Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id(
                             H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)),
                             "Unable to create nonlinear node memory type."),
                         H5Tclose};
  const auto insert = [&](const hid_t type, const hid_t integer_type,
                          const hid_t coordinate_type) {
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "internal_node_id", HOFFSET(NodeWriteRow, internal_node_id),
                  integer_type),
        "Unable to define nonlinear node ID.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "instance_name", HOFFSET(NodeWriteRow, instance_name),
                  string_type.Get()),
        "Unable to define nonlinear node instance.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "source_label", HOFFSET(NodeWriteRow, source_label),
                  string_type.Get()),
        "Unable to define nonlinear node source label.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "coordinates", HOFFSET(NodeWriteRow, coordinates),
                  coordinate_type),
        "Unable to define nonlinear node coordinates.");
  };
  insert(file_type.Get(), H5T_STD_U64LE, file_coordinates.Get());
  insert(memory_type.Get(), H5T_NATIVE_UINT64, memory_coordinates.Get());
  auto dataset = hdf5_internal::WriteCompoundDataset(
      file, "/model/nodes", rows.size(), file_type.Get(), memory_type.Get(),
      rows.data());
  hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U);
  hdf5_internal::WriteStringAttribute(dataset.Get(), "coordinate_system",
                                      "global-cartesian");
  hdf5_internal::WriteStringAttribute(dataset.Get(), "component_names", "X,Y,Z");
  hdf5_internal::WriteStringAttribute(dataset.Get(), "component_unit_dimensions",
                                      "length,length,length");
  hdf5_internal::WriteStringAttribute(dataset.Get(), "location",
                                      "nodal-reference");
}

void WriteElements(const hid_t file, const Domain& domain) {
  std::vector<ElementWriteRow> rows;
  rows.reserve(domain.ShellElements().Size());
  for (std::size_t index = 0U; index < domain.ShellElements().Size(); ++index) {
    const auto& element = domain.ShellElements()[index];
    rows.push_back({static_cast<std::uint64_t>(index),
                    element.source_id.instance_name.c_str(),
                    element.source_id.source_label_text.c_str(),
                    hdf5_internal::ShellSourceTypeName(element.source_type),
                    "FESA-MITC4-total-lagrangian",
                    {static_cast<std::uint64_t>(element.node_indices[0U]),
                     static_cast<std::uint64_t>(element.node_indices[1U]),
                     static_cast<std::uint64_t>(element.node_indices[2U]),
                     static_cast<std::uint64_t>(element.node_indices[3U])},
                    static_cast<std::uint64_t>(element.section_index),
                    static_cast<std::uint64_t>(element.material_index)});
  }
  auto string_type = hdf5_internal::MakeUtf8StringType();
  const hsize_t nodes[] = {4U};
  Hdf5Handle file_nodes{hdf5_internal::RequireHdf5Id(
                            H5Tarray_create2(H5T_STD_U64LE, 1, nodes),
                            "Unable to create shell connectivity type."),
                        H5Tclose};
  Hdf5Handle memory_nodes{hdf5_internal::RequireHdf5Id(
                              H5Tarray_create2(H5T_NATIVE_UINT64, 1, nodes),
                              "Unable to create shell connectivity memory type."),
                          H5Tclose};
  Hdf5Handle file_type{hdf5_internal::RequireHdf5Id(
                           H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)),
                           "Unable to create shell element file type."),
                       H5Tclose};
  Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id(
                             H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)),
                             "Unable to create shell element memory type."),
                         H5Tclose};
  const auto insert = [&](const hid_t type, const hid_t integer_type,
                          const hid_t node_type) {
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "internal_element_id",
                  HOFFSET(ElementWriteRow, internal_element_id), integer_type),
        "Unable to define shell element ID.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "instance_name", HOFFSET(ElementWriteRow, instance_name),
                  string_type.Get()),
        "Unable to define shell instance identity.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "source_label", HOFFSET(ElementWriteRow, source_label),
                  string_type.Get()),
        "Unable to define shell source label identity.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "source_element_type",
                  HOFFSET(ElementWriteRow, source_element_type),
                  string_type.Get()),
        "Unable to define shell source type identity.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "internal_formulation",
                  HOFFSET(ElementWriteRow, internal_formulation),
                  string_type.Get()),
        "Unable to define shell formulation identity.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "node_internal_ids", HOFFSET(ElementWriteRow, node_internal_ids),
                  node_type),
        "Unable to define shell connectivity.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "shell_section_internal_id",
                  HOFFSET(ElementWriteRow, shell_section_internal_id), integer_type),
        "Unable to define shell section identity.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "material_internal_id",
                  HOFFSET(ElementWriteRow, material_internal_id), integer_type),
        "Unable to define shell material identity.");
  };
  insert(file_type.Get(), H5T_STD_U64LE, file_nodes.Get());
  insert(memory_type.Get(), H5T_NATIVE_UINT64, memory_nodes.Get());
  auto dataset = hdf5_internal::WriteCompoundDataset(
      file, "/model/elements", rows.size(), file_type.Get(), memory_type.Get(),
      rows.data());
  hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U);
  hdf5_internal::WriteStringAttribute(dataset.Get(), "formulation",
                                      "FESA-MITC4-total-lagrangian");
}

void WriteMaterials(const hid_t file, const Domain& domain) {
  std::vector<MaterialWriteRow> rows;
  rows.reserve(domain.LinearElasticMaterials().Size());
  for (std::size_t index = 0U; index < domain.LinearElasticMaterials().Size();
       ++index) {
    const auto& material = domain.LinearElasticMaterials()[index];
    rows.push_back({static_cast<std::uint64_t>(index), material.name.c_str(),
                    material.youngs_modulus, material.poisson_ratio});
  }
  auto string_type = hdf5_internal::MakeUtf8StringType();
  Hdf5Handle file_type{hdf5_internal::RequireHdf5Id(
                           H5Tcreate(H5T_COMPOUND, sizeof(MaterialWriteRow)),
                           "Unable to create shell material file type."),
                       H5Tclose};
  Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id(
                             H5Tcreate(H5T_COMPOUND, sizeof(MaterialWriteRow)),
                             "Unable to create shell material memory type."),
                         H5Tclose};
  const auto insert = [&](const hid_t type, const hid_t integer_type,
                          const hid_t double_type) {
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "internal_material_id",
                  HOFFSET(MaterialWriteRow, internal_material_id), integer_type),
        "Unable to define shell material ID.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "name", HOFFSET(MaterialWriteRow, name), string_type.Get()),
        "Unable to define shell material name.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "E", HOFFSET(MaterialWriteRow, youngs_modulus),
                  double_type),
        "Unable to define shell material E.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "nu", HOFFSET(MaterialWriteRow, poisson_ratio),
                  double_type),
        "Unable to define shell material nu.");
  };
  insert(file_type.Get(), H5T_STD_U64LE, H5T_IEEE_F64LE);
  insert(memory_type.Get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE);
  auto dataset = hdf5_internal::WriteCompoundDataset(
      file, "/model/shell/materials", rows.size(), file_type.Get(),
      memory_type.Get(), rows.data());
  hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U);
  hdf5_internal::WriteStringAttribute(dataset.Get(), "component_unit_dimensions",
                                      "force/length^2,1");
  hdf5_internal::WriteStringAttribute(dataset.Get(), "location",
                                      "model-material");
}

void WriteSections(const hid_t file, const Domain& domain) {
  std::vector<std::string> source_files;
  source_files.reserve(domain.ShellSections().Size());
  for (std::size_t index = 0U; index < domain.ShellSections().Size(); ++index) {
    source_files.push_back(hdf5_internal::NormalizedPathString(
        domain.ShellSections()[index].location.file));
  }
  std::vector<SectionWriteRow> rows;
  rows.reserve(domain.ShellSections().Size());
  for (std::size_t index = 0U; index < domain.ShellSections().Size(); ++index) {
    const auto& section = domain.ShellSections()[index];
    rows.push_back({static_cast<std::uint64_t>(index), source_files[index].c_str(),
                    static_cast<std::uint64_t>(section.location.line),
                    section.name.c_str(),
                    static_cast<std::uint64_t>(section.material_index),
                    section.thickness});
  }
  auto string_type = hdf5_internal::MakeUtf8StringType();
  Hdf5Handle file_type{hdf5_internal::RequireHdf5Id(
                           H5Tcreate(H5T_COMPOUND, sizeof(SectionWriteRow)),
                           "Unable to create shell section file type."),
                       H5Tclose};
  Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id(
                             H5Tcreate(H5T_COMPOUND, sizeof(SectionWriteRow)),
                             "Unable to create shell section memory type."),
                         H5Tclose};
  const auto insert = [&](const hid_t type, const hid_t integer_type,
                          const hid_t double_type) {
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "internal_section_id",
                  HOFFSET(SectionWriteRow, internal_section_id), integer_type),
        "Unable to define shell section ID.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "source_file", HOFFSET(SectionWriteRow, source_file),
                  string_type.Get()),
        "Unable to define shell section source file.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "source_line", HOFFSET(SectionWriteRow, source_line),
                  integer_type),
        "Unable to define shell section source line.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "source_elset", HOFFSET(SectionWriteRow, source_elset),
                  string_type.Get()),
        "Unable to define shell section source elset.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "material_internal_id",
                  HOFFSET(SectionWriteRow, material_internal_id), integer_type),
        "Unable to define shell section material identity.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "thickness", HOFFSET(SectionWriteRow, thickness),
                  double_type),
        "Unable to define shell section thickness.");
  };
  insert(file_type.Get(), H5T_STD_U64LE, H5T_IEEE_F64LE);
  insert(memory_type.Get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE);
  auto dataset = hdf5_internal::WriteCompoundDataset(
      file, "/model/shell/sections", rows.size(), file_type.Get(),
      memory_type.Get(), rows.data());
  hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U);
  hdf5_internal::WriteStringAttribute(dataset.Get(), "component_unit_dimensions",
                                      "length");
  hdf5_internal::WriteStringAttribute(dataset.Get(), "location",
                                      "model-shell-section");
  hdf5_internal::WriteStringAttribute(dataset.Get(), "layering",
                                      "centered-single-layer");
}

void WriteModelData(const hid_t file, const Domain& domain,
                    const NonlinearModelData& model_data) {
  std::vector<double> directors;
  std::vector<double> frames;
  directors.reserve(domain.Nodes().size() * 3U);
  frames.reserve(domain.Nodes().size() * 9U);
  for (const auto& frame : domain.ShellNodeInitialFrames()) {
    directors.insert(directors.end(), frame.director.begin(), frame.director.end());
    frames.insert(frames.end(), frame.tangent_a.begin(), frame.tangent_a.end());
    frames.insert(frames.end(), frame.tangent_b.begin(), frame.tangent_b.end());
    frames.insert(frames.end(), frame.director.begin(), frame.director.end());
  }
  const hsize_t node_count = static_cast<hsize_t>(domain.Nodes().size());
  (void)WriteModelDoubleDataset(file, "/model/shell/nodal_director",
                                {node_count, 3U}, directors.data(),
                                directors.size(), "D1,D2,D3", "1,1,1",
                                "global-cartesian", "nodal-initial-director");
  auto initial_frame = WriteModelDoubleDataset(
      file, "/model/shell/nodal_frame", {node_count, 3U, 3U}, frames.data(),
      frames.size(), "X,Y,Z", "1,1,1", "global-cartesian",
      "nodal-initial-frame");
  hdf5_internal::WriteStringAttribute(initial_frame.Get(), "axis_names", "A,B,D");

  WriteMaterials(file, domain);
  WriteSections(file, domain);
  const std::vector<hsize_t> nodal_dimensions{node_count, kDofsPerNode};
  hdf5_internal::WriteUint8Dataset(file, "/model/nodal_constraint_mask",
                                   nodal_dimensions, model_data.constraint_mask.data(),
                                   model_data.constraint_mask.size());
  Hdf5Handle mask{hdf5_internal::RequireHdf5Id(
                      H5Dopen2(file, "/model/nodal_constraint_mask", H5P_DEFAULT),
                      "Unable to reopen nonlinear constraint mask."),
                  H5Dclose};
  hdf5_internal::WriteUint64Attribute(mask.Get(), "schema_version", 1U);
  hdf5_internal::WriteStringAttribute(mask.Get(), "component_unit_dimensions",
                                      "1,1,1,1,1,1");
  hdf5_internal::WriteStringAttribute(mask.Get(), "coordinate_system",
                                      "global-cartesian");
  hdf5_internal::WriteStringAttribute(mask.Get(), "location",
                                      "nodal-constraint");
  auto start = WriteModelDoubleDataset(
      file, "/model/prescribed_displacement_start", nodal_dimensions,
      model_data.prescribed_start.data(), model_data.prescribed_start.size(),
      "UX,UY,UZ,URX,URY,URZ", "length,length,length,1,1,1",
      "global-cartesian", "nodal-prescribed-start");
  hdf5_internal::WriteDoubleAttribute(start.Get(), "load_factor", 0.0);
  auto target = WriteModelDoubleDataset(
      file, "/model/prescribed_displacement_target", nodal_dimensions,
      model_data.prescribed_target.data(), model_data.prescribed_target.size(),
      "UX,UY,UZ,URX,URY,URZ", "length,length,length,1,1,1",
      "global-cartesian", "nodal-prescribed-target");
  hdf5_internal::WriteDoubleAttribute(target.Get(), "load_factor", 1.0);

  const double gauss = 1.0 / std::sqrt(3.0);
  const std::array<double, 8> locations{-gauss, -gauss, gauss,  -gauss,
                                        gauss,  gauss,  -gauss, gauss};
  (void)WriteModelDoubleDataset(file, "/model/shell/midsurface_locations",
                                {4U, 2U}, locations.data(), locations.size(),
                                "XI,ETA", "1,1", "shell-natural",
                                "midsurface-location");
  const std::array<double, 3> positions{-1.0, 0.0, 1.0};
  auto section_positions = WriteModelDoubleDataset(
      file, "/model/shell/section_positions", {3U}, positions.data(),
      positions.size(), "ZETA", "1", "shell-natural", "section-position");
  hdf5_internal::WriteStringAttribute(section_positions.Get(), "position_names",
                                      "BOTTOM,MIDDLE,TOP");
}

void WriteAcceptedFrameTable(const hid_t file,
                             const std::vector<AcceptedFrameState>& frames) {
  std::vector<AcceptedFrameWriteRow> rows;
  rows.reserve(frames.size());
  for (const auto& frame : frames) {
    rows.push_back({frame.frame_index, frame.increment_index, frame.load_factor,
                    frame.iteration_count,
                    frame.is_factor_one ? static_cast<std::uint8_t>(1U)
                                        : static_cast<std::uint8_t>(0U)});
  }
  Hdf5Handle file_type{hdf5_internal::RequireHdf5Id(
                           H5Tcreate(H5T_COMPOUND, sizeof(AcceptedFrameWriteRow)),
                           "Unable to create accepted-frame file type."),
                       H5Tclose};
  Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id(
                             H5Tcreate(H5T_COMPOUND, sizeof(AcceptedFrameWriteRow)),
                             "Unable to create accepted-frame memory type."),
                         H5Tclose};
  const auto insert = [&](const hid_t type, const hid_t integer_type,
                          const hid_t double_type, const hid_t flag_type) {
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "frame_index", HOFFSET(AcceptedFrameWriteRow, frame_index),
                  integer_type),
        "Unable to define accepted-frame index.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "increment_index",
                  HOFFSET(AcceptedFrameWriteRow, increment_index), integer_type),
        "Unable to define accepted increment index.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "load_factor", HOFFSET(AcceptedFrameWriteRow, load_factor),
                  double_type),
        "Unable to define accepted-frame load factor.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "iteration_count",
                  HOFFSET(AcceptedFrameWriteRow, iteration_count), integer_type),
        "Unable to define accepted-frame iteration count.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "is_factor_one",
                  HOFFSET(AcceptedFrameWriteRow, is_factor_one), flag_type),
        "Unable to define accepted-frame factor-one flag.");
  };
  insert(file_type.Get(), H5T_STD_U64LE, H5T_IEEE_F64LE, H5T_STD_U8LE);
  insert(memory_type.Get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE,
         H5T_NATIVE_UINT8);
  auto dataset = hdf5_internal::WriteCompoundDataset(
      file, "/steps/Step-1/accepted_frames", rows.size(), file_type.Get(),
      memory_type.Get(), rows.data());
  hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U);
}

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<double> NodalDisplacement(const AcceptedFrameState& frame) {
  std::vector<double> values(frame.displacement.Size());
  for (std::size_t node = 0U; node < frame.nodal_internal_order.size(); ++node) {
    const std::size_t offset = kDofsPerNode * node;
    for (std::size_t component = 0U; component < 3U; ++component) {
      values[offset + component] = frame.displacement[offset + component];
      values[offset + 3U + component] =
          frame.rotation_vector[offset + 3U + component];
    }
  }
  return values;
}

std::vector<double> FlattenTriads(const AcceptedFrameState& frame) {
  std::vector<double> values;
  values.reserve(frame.nodal_triad.size() * 9U);
  for (const auto& triad : frame.nodal_triad) {
    for (const auto& axis : triad) {
      values.insert(values.end(), axis.begin(), axis.end());
    }
  }
  return values;
}

std::vector<double> FlattenDirectors(const AcceptedFrameState& frame) {
  std::vector<double> values;
  values.reserve(frame.nodal_director.size() * 3U);
  for (const auto& director : frame.nodal_director) {
    values.insert(values.end(), director.begin(), director.end());
  }
  return values;
}

std::vector<double> SpatialReaction(const AcceptedFrameState& frame) {
  std::vector<double> values = CopyVector(frame.reaction);
  for (std::size_t node = 0U; node < frame.nodal_internal_order.size(); ++node) {
    const std::size_t offset = kDofsPerNode * node;
    const Vector3 phi{frame.rotation_vector[offset + 3U],
                      frame.rotation_vector[offset + 4U],
                      frame.rotation_vector[offset + 5U]};
    const Matrix transform = InverseTransposeLeftJacobianSO3(phi);
    const std::array<double, 3> generalized{
        frame.reaction[offset + 3U], frame.reaction[offset + 4U],
        frame.reaction[offset + 5U]};
    for (std::size_t row = 0U; row < 3U; ++row) {
      values[offset + 3U + row] = 0.0;
      for (std::size_t column = 0U; column < 3U; ++column) {
        values[offset + 3U + row] += transform(row, column) * generalized[column];
      }
    }
  }
  return values;
}

void WriteShellIdentity(const hid_t dataset, const bool has_section_positions,
                        const std::string& local_frame_path) {
  hdf5_internal::WriteStringAttribute(dataset, "source_element_type_dataset",
                                      "/model/elements.source_element_type");
  hdf5_internal::WriteStringAttribute(dataset, "internal_formulation",
                                      "FESA-MITC4-total-lagrangian");
  hdf5_internal::WriteStringAttribute(dataset, "midsurface_location_dataset",
                                      "/model/shell/midsurface_locations");
  hdf5_internal::WriteStringAttribute(dataset, "local_frame_dataset",
                                      local_frame_path);
  if (has_section_positions) {
    hdf5_internal::WriteStringAttribute(dataset, "section_position_dataset",
                                        "/model/shell/section_positions");
  }
}

void WriteFrame(const hid_t file, const Domain& domain,
                const AcceptedFrameState& frame) {
  const std::string root = "/steps/Step-1/frames/" +
                           std::to_string(frame.frame_index);
  (void)hdf5_internal::CreateGroup(file, (root + "/nodal").c_str());
  (void)hdf5_internal::CreateGroup(file, (root + "/element/shell").c_str());
  (void)hdf5_internal::CreateGroup(file, (root + "/global").c_str());
  const hsize_t node_count = static_cast<hsize_t>(domain.Nodes().size());
  const std::vector<hsize_t> nodal_dimensions{node_count, kDofsPerNode};
  const auto displacement = NodalDisplacement(frame);
  (void)WriteFrameDoubleDataset(
      file, root + "/nodal/displacement", nodal_dimensions, displacement.data(),
      displacement.size(), "UX,UY,UZ,URX,URY,URZ", "length,length,length,1,1,1",
      "global-cartesian", "nodal", frame);
  const auto triads = FlattenTriads(frame);
  (void)WriteFrameDoubleDataset(file, root + "/nodal/triad",
                                {node_count, 3U, 3U}, triads.data(), triads.size(),
                                "X,Y,Z", "1,1,1", "global-cartesian",
                                "nodal-triad", frame);
  const auto directors = FlattenDirectors(frame);
  (void)WriteFrameDoubleDataset(file, root + "/nodal/director", {node_count, 3U},
                                directors.data(), directors.size(), "D1,D2,D3",
                                "1,1,1", "global-cartesian", "nodal-director",
                                frame);
  const auto internal = CopyVector(frame.internal_force);
  (void)WriteFrameDoubleDataset(
      file, root + "/nodal/generalized_internal_force", nodal_dimensions,
      internal.data(), internal.size(), "FI1,FI2,FI3,MI_PHI1,MI_PHI2,MI_PHI3",
      "force,force,force,force*length,force*length,force*length",
      "global-cartesian", "nodal", frame);
  const auto external = CopyVector(frame.external_force);
  (void)WriteFrameDoubleDataset(
      file, root + "/nodal/generalized_external_force", nodal_dimensions,
      external.data(), external.size(), "FE1,FE2,FE3,ME_PHI1,ME_PHI2,ME_PHI3",
      "force,force,force,force*length,force*length,force*length",
      "global-cartesian", "nodal", frame);
  const auto residual = CopyVector(frame.residual);
  (void)WriteFrameDoubleDataset(
      file, root + "/nodal/generalized_residual", nodal_dimensions,
      residual.data(), residual.size(), "RU1,RU2,RU3,RPHI1,RPHI2,RPHI3",
      "force,force,force,force*length,force*length,force*length",
      "global-cartesian", "nodal", frame);
  const auto reaction = SpatialReaction(frame);
  (void)WriteFrameDoubleDataset(
      file, root + "/nodal/reaction", nodal_dimensions, reaction.data(),
      reaction.size(), "RF1,RF2,RF3,RM1,RM2,RM3",
      "force,force,force,force*length,force*length,force*length",
      "global-cartesian", "nodal", frame);

  std::vector<double> local_frames;
  std::vector<double> strains;
  std::vector<double> resultants;
  std::vector<double> stresses;
  local_frames.reserve(frame.shell_results.size() * 9U);
  strains.reserve(frame.shell_results.size() * 8U);
  resultants.reserve(frame.shell_results.size() * 8U);
  stresses.reserve(frame.shell_results.size() * kSectionPositions * 3U);
  for (const auto& row : frame.shell_results) {
    for (const auto& axis : row.local_frame) {
      local_frames.insert(local_frames.end(), axis.begin(), axis.end());
    }
    strains.insert(strains.end(), row.generalized_strain.begin(),
                   row.generalized_strain.end());
    resultants.insert(resultants.end(), row.section_resultant.begin(),
                      row.section_resultant.end());
    for (const auto& stress : row.stress) {
      stresses.insert(stresses.end(), stress.components.begin(),
                      stress.components.end());
    }
  }
  const hsize_t element_count = static_cast<hsize_t>(domain.ShellElements().Size());
  const std::string shell_root = root + "/element/shell";
  auto local_frame = WriteFrameDoubleDataset(
      file, shell_root + "/local_frame", {element_count, 4U, 3U, 3U},
      local_frames.data(), local_frames.size(), "X,Y,Z", "1,1,1",
      "global-cartesian", "shell-local-frame", frame);
  hdf5_internal::WriteStringAttribute(local_frame.Get(), "axis_names", "E1,E2,E3");
  hdf5_internal::WriteStringAttribute(local_frame.Get(), "source_element_type_dataset",
                                      "/model/elements.source_element_type");
  hdf5_internal::WriteStringAttribute(local_frame.Get(), "internal_formulation",
                                      "FESA-MITC4-total-lagrangian");
  hdf5_internal::WriteStringAttribute(local_frame.Get(), "midsurface_location_dataset",
                                      "/model/shell/midsurface_locations");
  auto strain = WriteFrameDoubleDataset(
      file, shell_root + "/generalized_strain", {element_count, 4U, 8U},
      strains.data(), strains.size(), "E11,E22,G12,K11,K22,K12,G13,G23",
      "1,1,1,1/length,1/length,1/length,1,1", "shell-local", "midsurface",
      frame);
  WriteShellIdentity(strain.Get(), false, shell_root + "/local_frame");
  auto resultant = WriteFrameDoubleDataset(
      file, shell_root + "/section_resultant", {element_count, 4U, 8U},
      resultants.data(), 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", frame);
  WriteShellIdentity(resultant.Get(), false, shell_root + "/local_frame");
  auto stress = WriteFrameDoubleDataset(
      file, shell_root + "/stress", {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", frame);
  WriteShellIdentity(stress.Get(), true, shell_root + "/local_frame");
  (void)WriteFrameDoubleDataset(file, root + "/global/energy", {1U},
                                &frame.physical_strain_energy, 1U,
                                "PHYSICAL_STRAIN_ENERGY", "force*length",
                                "global", "global", frame);
  (void)WriteFrameDoubleDataset(
      file, root + "/global/equilibrium", {6U}, frame.global_equilibrium.data(),
      frame.global_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-balance-error", frame);
}

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(hdf5_internal::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 = hdf5_internal::MakeUtf8StringType();
  Hdf5Handle file_type{hdf5_internal::RequireHdf5Id(
                           H5Tcreate(H5T_COMPOUND, sizeof(DiagnosticWriteRow)),
                           "Unable to create nonlinear diagnostic file type."),
                       H5Tclose};
  Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id(
                             H5Tcreate(H5T_COMPOUND, sizeof(DiagnosticWriteRow)),
                             "Unable to create nonlinear diagnostic memory type."),
                         H5Tclose};
  const auto insert = [&](const hid_t type, const hid_t integer_type) {
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "severity", HOFFSET(DiagnosticWriteRow, severity),
                  string_type.Get()),
        "Unable to define nonlinear diagnostic severity.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "code", HOFFSET(DiagnosticWriteRow, code),
                  string_type.Get()),
        "Unable to define nonlinear diagnostic code.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "file", HOFFSET(DiagnosticWriteRow, file),
                  string_type.Get()),
        "Unable to define nonlinear diagnostic file.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "line", HOFFSET(DiagnosticWriteRow, line), integer_type),
        "Unable to define nonlinear diagnostic line.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "keyword", HOFFSET(DiagnosticWriteRow, keyword),
                  string_type.Get()),
        "Unable to define nonlinear diagnostic keyword.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "entity_identity",
                  HOFFSET(DiagnosticWriteRow, entity_identity), string_type.Get()),
        "Unable to define nonlinear diagnostic entity identity.");
    hdf5_internal::RequireHdf5(
        H5Tinsert(type, "message", HOFFSET(DiagnosticWriteRow, message),
                  string_type.Get()),
        "Unable to define nonlinear diagnostic message.");
  };
  insert(file_type.Get(), H5T_STD_U64LE);
  insert(memory_type.Get(), H5T_NATIVE_UINT64);
  auto dataset = hdf5_internal::WriteCompoundDataset(
      file, "/diagnostics", rows.size(), file_type.Get(), memory_type.Get(),
      rows.data());
  hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U);
}

Status ValidateDiagnostics(const std::vector<Diagnostic>& diagnostics) {
  for (const auto& diagnostic : diagnostics) {
    if (!hdf5_internal::IsValidUtf8(diagnostic.code) ||
        !hdf5_internal::IsValidUtf8(diagnostic.keyword) ||
        !hdf5_internal::IsValidUtf8(diagnostic.entity_identity) ||
        !hdf5_internal::IsValidUtf8(diagnostic.message)) {
      return Failure("invalid-result-diagnostic",
                     "Nonlinear diagnostic text must be valid UTF-8.");
    }
  }
  return Status::Ok();
}

}  // namespace

Status NonlinearHdf5ResultsWriter::Write(
    const std::filesystem::path& output_path, const Domain& domain,
    const NonlinearAnalysisState& state,
    const std::vector<Diagnostic>& diagnostics) {
  if (output_path.empty() || output_path.filename().empty()) {
    return Failure("invalid-output-path", "The HDF5 output path must name a file.");
  }
  const Status domain_status = ValidateDomain(domain);
  if (!domain_status.IsOk()) {
    return domain_status;
  }
  const Status state_status = ValidateState(domain, state);
  if (!state_status.IsOk()) {
    return state_status;
  }
  const Status diagnostic_status = ValidateDiagnostics(diagnostics);
  if (!diagnostic_status.IsOk()) {
    return diagnostic_status;
  }
  NonlinearModelData model_data;
  const Status model_status = BuildModelData(domain, model_data);
  if (!model_status.IsOk()) {
    return model_status;
  }

  try {
    hdf5_internal::Hdf5ErrorSilencer silence_backend_errors;
    const std::filesystem::path candidate_path =
        hdf5_internal::MakeCandidatePath(output_path);
    hdf5_internal::CandidateFileGuard cleanup{candidate_path};
    Hdf5Handle file{hdf5_internal::RequireHdf5Id(
                        H5Fcreate(candidate_path.string().c_str(), H5F_ACC_EXCL,
                                  H5P_DEFAULT, H5P_DEFAULT),
                        "Unable to create nonlinear temporary HDF5 file."),
                    H5Fclose};
    (void)hdf5_internal::CreateGroup(file.Get(), "/model");
    (void)hdf5_internal::CreateGroup(file.Get(), "/model/shell");
    (void)hdf5_internal::CreateGroup(file.Get(), "/steps/Step-1");
    WriteMetadata(file.Get(), domain);
    WriteNodes(file.Get(), domain);
    WriteElements(file.Get(), domain);
    WriteModelData(file.Get(), domain, model_data);
    WriteAcceptedFrameTable(file.Get(), state.AcceptedFrames());
    for (const auto& frame : state.AcceptedFrames()) {
      WriteFrame(file.Get(), domain, frame);
    }
    WriteDiagnostics(file.Get(), diagnostics);
    hdf5_internal::RequireHdf5(H5Fflush(file.Get(), H5F_SCOPE_GLOBAL),
                               "Unable to flush nonlinear temporary HDF5 file.");
    hdf5_internal::RequireHdf5(
        file.CloseChecked(), "Unable to close nonlinear temporary HDF5 file.");
    const Status self_check = hdf5_internal::NonlinearSelfCheck(
        candidate_path, domain, state, diagnostics.size());
    const Status atomic = hdf5_internal::AtomicFinalizeValidatedCandidate(
        candidate_path, output_path, self_check);
    if (atomic.IsOk()) {
      cleanup.Release();
    }
    return atomic;
  } catch (const Hdf5Failure& failure) {
    return Failure("hdf5-write-failure", failure.what());
  } catch (const std::exception& failure) {
    return Failure("hdf5-write-failure", failure.what());
  }
}

}  // namespace fesa