feat(cpp-object-oriented-modular-refactoring): step 23 - hdf5-writer-modules

src/fesa/io/hdf5/hdf5_result_writer.cpp

@@ -0,0 +1,955 @@
+#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