diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index cbb43b6..6013f2b 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -25,7 +25,12 @@ add_library( io/abaqus/material_property_mapper.cpp io/abaqus/step_definition_mapper.cpp io/abaqus/topology_mapper.cpp + io/hdf5/hdf5_atomic_file.cpp + io/hdf5/hdf5_model_writer.cpp + io/hdf5/hdf5_primitives.cpp + io/hdf5/hdf5_result_writer.cpp io/hdf5/hdf5_results_writer.cpp + io/hdf5/hdf5_self_check.cpp loads/concentrated_nodal_load.cpp materials/isotropic_linear_elastic_material.cpp math/dense_blas_internal.cpp diff --git a/src/fesa/io/hdf5/hdf5_atomic_file.cpp b/src/fesa/io/hdf5/hdf5_atomic_file.cpp new file mode 100644 index 0000000..4bf91f5 --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_atomic_file.cpp @@ -0,0 +1,88 @@ +#define NOMINMAX +#include "io/hdf5/hdf5_atomic_file.h" + +#include + +#include +#include +#include +#include +#include +#include + +#include "io/hdf5/hdf5_primitives.h" +#include "io/hdf5/hdf5_raii.h" + +namespace fesa::hdf5_internal { +namespace { + +bool AtomicReplace(const std::filesystem::path& candidate_path, + const std::filesystem::path& output_path) { + const DWORD attributes = GetFileAttributesW(output_path.c_str()); + if (attributes != INVALID_FILE_ATTRIBUTES) { + return ReplaceFileW(output_path.c_str(), candidate_path.c_str(), nullptr, + 0U, nullptr, nullptr) != FALSE; + } + const DWORD error = GetLastError(); + if (error != ERROR_FILE_NOT_FOUND && error != ERROR_PATH_NOT_FOUND) { + return false; + } + return MoveFileExW(candidate_path.c_str(), output_path.c_str(), + MOVEFILE_WRITE_THROUGH) != FALSE; +} + +} // namespace + +CandidateFileGuard::CandidateFileGuard(std::filesystem::path path) + : path_{std::move(path)} {} + +CandidateFileGuard::~CandidateFileGuard() { + if (active_) { + std::error_code ignored; + (void)std::filesystem::remove(path_, ignored); + } +} + +void CandidateFileGuard::Release() noexcept { active_ = false; } + +std::filesystem::path MakeCandidatePath( + const std::filesystem::path& output_path) { + static std::atomic sequence{0U}; + const std::filesystem::path parent = output_path.parent_path(); + for (std::size_t attempt = 0U; attempt < 1024U; ++attempt) { + const std::wstring filename = + L"." + output_path.filename().wstring() + L".tmp." + + std::to_wstring(GetCurrentProcessId()) + L"." + + std::to_wstring(sequence.fetch_add(1U)); + const std::filesystem::path candidate = parent / filename; + std::error_code error; + const bool exists = std::filesystem::exists(candidate, error); + if (error) { + throw Hdf5Failure{"Unable to inspect the temporary output path."}; + } + if (!exists) { + return candidate; + } + } + throw Hdf5Failure{"Unable to allocate a unique temporary output path."}; +} + +Status AtomicFinalizeValidatedCandidate( + const std::filesystem::path& candidate_path, + const std::filesystem::path& output_path, const Status& validation) { + if (!validation.IsOk()) { + std::error_code ignored; + (void)std::filesystem::remove(candidate_path, ignored); + return validation; + } + if (!AtomicReplace(candidate_path, output_path)) { + std::error_code ignored; + (void)std::filesystem::remove(candidate_path, ignored); + return OutputFailure("hdf5-finalization-failure", + "The checked temporary HDF5 file could not replace " + "the final output."); + } + return Status::Ok(); +} + +} // namespace fesa::hdf5_internal diff --git a/src/fesa/io/hdf5/hdf5_atomic_file.h b/src/fesa/io/hdf5/hdf5_atomic_file.h new file mode 100644 index 0000000..98d540d --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_atomic_file.h @@ -0,0 +1,37 @@ +#ifndef FESA_IO_HDF5_HDF5_ATOMIC_FILE_H_ +#define FESA_IO_HDF5_HDF5_ATOMIC_FILE_H_ + +#include + +#include "fesa/core/status.h" + +namespace fesa::hdf5_internal { + +/// @brief Removes an uncommitted same-directory candidate on failure. +class CandidateFileGuard { + public: + explicit CandidateFileGuard(std::filesystem::path path); + CandidateFileGuard(const CandidateFileGuard&) = delete; + CandidateFileGuard& operator=(const CandidateFileGuard&) = delete; + ~CandidateFileGuard(); + + void Release() noexcept; + + private: + std::filesystem::path path_; + bool active_{true}; +}; + +/// @brief Selects a unique same-directory candidate path. +std::filesystem::path MakeCandidatePath( + const std::filesystem::path& output_path); + +/// @brief Publishes only a successfully self-checked candidate. +/// @note Validation failure removes the candidate and preserves any final. +Status AtomicFinalizeValidatedCandidate( + const std::filesystem::path& candidate_path, + const std::filesystem::path& output_path, const Status& validation); + +} // namespace fesa::hdf5_internal + +#endif // FESA_IO_HDF5_HDF5_ATOMIC_FILE_H_ diff --git a/src/fesa/io/hdf5/hdf5_model_writer.cpp b/src/fesa/io/hdf5/hdf5_model_writer.cpp new file mode 100644 index 0000000..628c331 --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_model_writer.cpp @@ -0,0 +1,491 @@ +#include "io/hdf5/hdf5_model_writer.h" + +#include +#include +#include +#include +#include + +#include "fesa/build_info.h" + +namespace fesa::hdf5_internal { + +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; + std::uint64_t node_internal_ids[2]; + double local_axes[9]; +}; + +struct ShellElementWriteRow { + 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 ShellMaterialWriteRow { + std::uint64_t internal_material_id; + const char* name; + double youngs_modulus; + double poisson_ratio; +}; + +struct ShellSectionWriteRow { + 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; +}; + +void WriteMetadata(const hid_t file, const Domain& domain) { + auto metadata = CreateGroup(file, "/metadata"); + WriteUint64Attribute(metadata.Get(), "schema_version", 0U); + WriteStringAttribute(metadata.Get(), "solver_version", + std::string{SolverVersion()}); + WriteStringAttribute(metadata.Get(), "source_input_identity", + SourceInputIdentity(domain)); + WriteStringAttribute(metadata.Get(), "unit_system_label", + "user-consistent-unspecified"); + if (IsShellDomain(domain)) { + WriteStringAttribute(metadata.Get(), "feature_id", + "linear-static-mitc4-shell"); + WriteStringAttribute( + metadata.Get(), "coordinate_convention", + "global-cartesian; shell-local=(e1,e2,e3); positive-thickness=+zeta"); + WriteStringAttribute(metadata.Get(), "internal_formulation", "FESA-MITC4"); + WriteStringAttribute(metadata.Get(), "integration_rule", + "2x2x2-gauss; mitc4-edge-midpoint-shear"); + } else { + WriteStringAttribute(metadata.Get(), "feature_id", + "linear-static-3d-euler-beam"); + WriteStringAttribute(metadata.Get(), "coordinate_convention", + "global-cartesian; beam-local=(t,n1,t-cross-n1)"); + WriteStringAttribute(metadata.Get(), "element_formulation", + "B33-3D-Euler-Bernoulli"); + } + WriteStringAttribute(metadata.Get(), "step_name", kStepName); + WriteUint64Attribute(metadata.Get(), "frame_index", 0U); +} + +void WriteNodes(const hid_t file, const Domain& domain) { + std::vector 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(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 = MakeUtf8StringType(); + const hsize_t coordinate_dimensions[] = {3U}; + Hdf5Handle file_coordinates{ + RequireHdf5Id(H5Tarray_create2(H5T_IEEE_F64LE, 1, coordinate_dimensions), + "Unable to create the node coordinate file type."), + H5Tclose}; + Hdf5Handle memory_coordinates{ + RequireHdf5Id( + H5Tarray_create2(H5T_NATIVE_DOUBLE, 1, coordinate_dimensions), + "Unable to create the node coordinate memory type."), + H5Tclose}; + Hdf5Handle file_type{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)), + "Unable to create the node file type."), + H5Tclose}; + Hdf5Handle memory_type{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)), + "Unable to create the node memory type."), + H5Tclose}; + RequireHdf5(H5Tinsert(file_type.Get(), "internal_node_id", + HOFFSET(NodeWriteRow, internal_node_id), H5T_STD_U64LE), + "Unable to define the node internal ID field."); + RequireHdf5( + H5Tinsert(file_type.Get(), "instance_name", + HOFFSET(NodeWriteRow, instance_name), string_type.Get()), + "Unable to define the node instance field."); + RequireHdf5(H5Tinsert(file_type.Get(), "source_label", + HOFFSET(NodeWriteRow, source_label), string_type.Get()), + "Unable to define the node source-label field."); + RequireHdf5( + H5Tinsert(file_type.Get(), "coordinates", + HOFFSET(NodeWriteRow, coordinates), file_coordinates.Get()), + "Unable to define the node coordinate field."); + RequireHdf5( + H5Tinsert(memory_type.Get(), "internal_node_id", + HOFFSET(NodeWriteRow, internal_node_id), H5T_NATIVE_UINT64), + "Unable to define the node internal ID memory field."); + RequireHdf5( + H5Tinsert(memory_type.Get(), "instance_name", + HOFFSET(NodeWriteRow, instance_name), string_type.Get()), + "Unable to define the node instance memory field."); + RequireHdf5(H5Tinsert(memory_type.Get(), "source_label", + HOFFSET(NodeWriteRow, source_label), string_type.Get()), + "Unable to define the node source-label memory field."); + RequireHdf5( + H5Tinsert(memory_type.Get(), "coordinates", + HOFFSET(NodeWriteRow, coordinates), memory_coordinates.Get()), + "Unable to define the node coordinate memory field."); + + auto dataset = + WriteCompoundDataset(file, "/model/nodes", rows.size(), file_type.Get(), + memory_type.Get(), rows.data()); + WriteStringAttribute(dataset.Get(), "coordinate_system", "global-cartesian"); + WriteStringAttribute(dataset.Get(), "units_label", "length"); +} + +void WriteBeamElements(const hid_t file, const Domain& domain, + const std::vector& axes) { + std::vector rows; + rows.reserve(domain.BeamElements().Size()); + for (std::size_t index = 0U; index < domain.BeamElements().Size(); ++index) { + const auto& element = domain.BeamElements()[index]; + ElementWriteRow row{static_cast(index), + element.source_id.instance_name.c_str(), + element.source_id.source_label_text.c_str(), + {static_cast(element.node_indices[0U]), + static_cast(element.node_indices[1U])}, + {}}; + std::copy(axes[index].begin(), axes[index].end(), row.local_axes); + rows.push_back(row); + } + + auto string_type = MakeUtf8StringType(); + const hsize_t node_dimensions[] = {2U}; + const hsize_t axis_dimensions[] = {3U, 3U}; + Hdf5Handle file_nodes{ + RequireHdf5Id(H5Tarray_create2(H5T_STD_U64LE, 1, node_dimensions), + "Unable to create the element connectivity file type."), + H5Tclose}; + Hdf5Handle memory_nodes{ + RequireHdf5Id(H5Tarray_create2(H5T_NATIVE_UINT64, 1, node_dimensions), + "Unable to create the element connectivity memory type."), + H5Tclose}; + Hdf5Handle file_axes{ + RequireHdf5Id(H5Tarray_create2(H5T_IEEE_F64LE, 2, axis_dimensions), + "Unable to create the local-axis file type."), + H5Tclose}; + Hdf5Handle memory_axes{ + RequireHdf5Id(H5Tarray_create2(H5T_NATIVE_DOUBLE, 2, axis_dimensions), + "Unable to create the local-axis memory type."), + H5Tclose}; + Hdf5Handle file_type{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)), + "Unable to create the element file type."), + H5Tclose}; + Hdf5Handle memory_type{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)), + "Unable to create the element memory type."), + H5Tclose}; + const auto insert_fields = [&](const hid_t type, const hid_t integer_type, + const hid_t node_type, const hid_t axes_type) { + RequireHdf5( + H5Tinsert(type, "internal_element_id", + HOFFSET(ElementWriteRow, internal_element_id), integer_type), + "Unable to define the element internal ID field."); + RequireHdf5( + H5Tinsert(type, "instance_name", + HOFFSET(ElementWriteRow, instance_name), string_type.Get()), + "Unable to define the element instance field."); + RequireHdf5( + H5Tinsert(type, "source_label", HOFFSET(ElementWriteRow, source_label), + string_type.Get()), + "Unable to define the element source-label field."); + RequireHdf5( + H5Tinsert(type, "node_internal_ids", + HOFFSET(ElementWriteRow, node_internal_ids), node_type), + "Unable to define the element connectivity field."); + RequireHdf5(H5Tinsert(type, "local_axes", + HOFFSET(ElementWriteRow, local_axes), axes_type), + "Unable to define the element local-axis field."); + }; + insert_fields(file_type.Get(), H5T_STD_U64LE, file_nodes.Get(), + file_axes.Get()); + insert_fields(memory_type.Get(), H5T_NATIVE_UINT64, memory_nodes.Get(), + memory_axes.Get()); + auto dataset = + WriteCompoundDataset(file, "/model/elements", rows.size(), + file_type.Get(), memory_type.Get(), rows.data()); + WriteStringAttribute(dataset.Get(), "formulation", "B33-3D-Euler-Bernoulli"); +} + +void WriteShellElements(const hid_t file, const Domain& domain) { + std::vector 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(index), + element.source_id.instance_name.c_str(), + element.source_id.source_label_text.c_str(), + ShellSourceTypeName(element.source_type), + kMitc4InternalFormulation.data(), + {static_cast(element.node_indices[0U]), + static_cast(element.node_indices[1U]), + static_cast(element.node_indices[2U]), + static_cast(element.node_indices[3U])}, + static_cast(element.section_index), + static_cast(element.material_index)}); + } + + auto string_type = MakeUtf8StringType(); + const hsize_t node_dimensions[] = {kShellNodeCount}; + Hdf5Handle file_nodes{ + RequireHdf5Id(H5Tarray_create2(H5T_STD_U64LE, 1, node_dimensions), + "Unable to create shell connectivity file type."), + H5Tclose}; + Hdf5Handle memory_nodes{ + RequireHdf5Id(H5Tarray_create2(H5T_NATIVE_UINT64, 1, node_dimensions), + "Unable to create shell connectivity memory type."), + H5Tclose}; + Hdf5Handle file_type{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellElementWriteRow)), + "Unable to create shell element file type."), + H5Tclose}; + Hdf5Handle memory_type{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellElementWriteRow)), + "Unable to create shell element memory type."), + H5Tclose}; + const auto insert_fields = [&](const hid_t type, const hid_t integer_type, + const hid_t node_type) { + RequireHdf5(H5Tinsert(type, "internal_element_id", + HOFFSET(ShellElementWriteRow, internal_element_id), + integer_type), + "Unable to define shell element ID field."); + RequireHdf5(H5Tinsert(type, "instance_name", + HOFFSET(ShellElementWriteRow, instance_name), + string_type.Get()), + "Unable to define shell element instance field."); + RequireHdf5(H5Tinsert(type, "source_label", + HOFFSET(ShellElementWriteRow, source_label), + string_type.Get()), + "Unable to define shell element source-label field."); + RequireHdf5(H5Tinsert(type, "source_element_type", + HOFFSET(ShellElementWriteRow, source_element_type), + string_type.Get()), + "Unable to define shell source-type field."); + RequireHdf5(H5Tinsert(type, "internal_formulation", + HOFFSET(ShellElementWriteRow, internal_formulation), + string_type.Get()), + "Unable to define shell formulation field."); + RequireHdf5( + H5Tinsert(type, "node_internal_ids", + HOFFSET(ShellElementWriteRow, node_internal_ids), node_type), + "Unable to define shell connectivity field."); + RequireHdf5( + H5Tinsert(type, "shell_section_internal_id", + HOFFSET(ShellElementWriteRow, shell_section_internal_id), + integer_type), + "Unable to define shell section ID field."); + RequireHdf5(H5Tinsert(type, "material_internal_id", + HOFFSET(ShellElementWriteRow, material_internal_id), + integer_type), + "Unable to define shell material ID field."); + }; + insert_fields(file_type.Get(), H5T_STD_U64LE, file_nodes.Get()); + insert_fields(memory_type.Get(), H5T_NATIVE_UINT64, memory_nodes.Get()); + auto dataset = + WriteCompoundDataset(file, "/model/elements", rows.size(), + file_type.Get(), memory_type.Get(), rows.data()); + WriteStringAttribute(dataset.Get(), "formulation", "FESA-MITC4"); +} + +void WriteShellMaterials(const hid_t file, const Domain& domain) { + std::vector 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(index), material.name.c_str(), + material.youngs_modulus, material.poisson_ratio}); + } + auto string_type = MakeUtf8StringType(); + Hdf5Handle file_type{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellMaterialWriteRow)), + "Unable to create shell material file type."), + H5Tclose}; + Hdf5Handle memory_type{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellMaterialWriteRow)), + "Unable to create shell material 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_material_id", + HOFFSET(ShellMaterialWriteRow, internal_material_id), + integer_type), + "Unable to define shell material ID field."); + RequireHdf5(H5Tinsert(type, "name", HOFFSET(ShellMaterialWriteRow, name), + string_type.Get()), + "Unable to define shell material name field."); + RequireHdf5( + H5Tinsert(type, "E", HOFFSET(ShellMaterialWriteRow, youngs_modulus), + double_type), + "Unable to define shell material E field."); + RequireHdf5( + H5Tinsert(type, "nu", HOFFSET(ShellMaterialWriteRow, poisson_ratio), + double_type), + "Unable to define shell material nu 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, "/model/shell/materials", rows.size(), + file_type.Get(), memory_type.Get(), rows.data()); + WriteStringAttribute(dataset.Get(), "component_unit_dimensions", + "force/length^2,1"); +} + +void WriteShellSections(const hid_t file, const Domain& domain) { + std::vector source_files; + source_files.reserve(domain.ShellSections().Size()); + for (std::size_t index = 0U; index < domain.ShellSections().Size(); ++index) { + const auto& section = domain.ShellSections()[index]; + source_files.push_back(NormalizedPathString(section.location.file)); + } + std::vector 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(index), + source_files[index].c_str(), + static_cast(section.location.line), + section.name.c_str(), + static_cast(section.material_index), + section.thickness}); + } + auto string_type = MakeUtf8StringType(); + Hdf5Handle file_type{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellSectionWriteRow)), + "Unable to create shell section file type."), + H5Tclose}; + Hdf5Handle memory_type{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellSectionWriteRow)), + "Unable to create shell section 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_section_id", + HOFFSET(ShellSectionWriteRow, internal_section_id), + integer_type), + "Unable to define shell section ID field."); + RequireHdf5(H5Tinsert(type, "source_file", + HOFFSET(ShellSectionWriteRow, source_file), + string_type.Get()), + "Unable to define shell section source-file field."); + RequireHdf5( + H5Tinsert(type, "source_line", + HOFFSET(ShellSectionWriteRow, source_line), integer_type), + "Unable to define shell section source-line field."); + RequireHdf5(H5Tinsert(type, "source_elset", + HOFFSET(ShellSectionWriteRow, source_elset), + string_type.Get()), + "Unable to define shell section ELSET field."); + RequireHdf5(H5Tinsert(type, "material_internal_id", + HOFFSET(ShellSectionWriteRow, material_internal_id), + integer_type), + "Unable to define shell section material ID field."); + RequireHdf5( + H5Tinsert(type, "thickness", HOFFSET(ShellSectionWriteRow, thickness), + double_type), + "Unable to define shell section thickness 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, "/model/shell/sections", rows.size(), + file_type.Get(), memory_type.Get(), rows.data()); + WriteStringAttribute(dataset.Get(), "thickness_unit_dimension", "length"); + WriteStringAttribute(dataset.Get(), "layering", "centered-single-layer"); +} + +void WriteShellModelData(const hid_t file, const Domain& domain, + const WriterModelData& model_data) { + std::vector directors; + directors.reserve(domain.Nodes().size() * 3U); + std::vector frames; + 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()); + } + WriteModelDoubleDataset(file, "/model/shell/nodal_director", + {static_cast(domain.Nodes().size()), 3U}, + directors.data(), directors.size(), "D1,D2,D3", + "1,1,1", "global-cartesian", "nodal"); + WriteModelDoubleDataset(file, "/model/shell/nodal_frame", + {static_cast(domain.Nodes().size()), 3U, 3U}, + frames.data(), frames.size(), "X,Y,Z", "1,1,1", + "global-cartesian", "nodal-frame"); + { + Hdf5Handle dataset{ + RequireHdf5Id(H5Dopen2(file, "/model/shell/nodal_frame", H5P_DEFAULT), + "Unable to reopen the nodal-frame dataset."), + H5Dclose}; + WriteStringAttribute(dataset.Get(), "axis_names", "A,B,D"); + } + + WriteShellMaterials(file, domain); + WriteShellSections(file, domain); + const std::vector nodal_dimensions{ + static_cast(domain.Nodes().size()), kDofsPerNode}; + WriteUint8Dataset(file, "/model/nodal_constraint_mask", nodal_dimensions, + model_data.constraint_mask.data(), + model_data.constraint_mask.size()); + WriteModelDoubleDataset( + file, "/model/prescribed_displacement", nodal_dimensions, + model_data.prescribed_displacement.data(), + model_data.prescribed_displacement.size(), "UX,UY,UZ,URX,URY,URZ", + "length,length,length,radian,radian,radian", "global-cartesian", + "nodal-prescribed-value"); + + const double gauss = 1.0 / std::sqrt(3.0); + const std::array locations{-gauss, -gauss, gauss, -gauss, + gauss, gauss, -gauss, gauss}; + WriteModelDoubleDataset(file, "/model/shell/midsurface_locations", {4U, 2U}, + locations.data(), locations.size(), "XI,ETA", "1,1", + "shell-natural", "midsurface-location"); + const std::array section_positions{-1.0, 0.0, 1.0}; + WriteModelDoubleDataset(file, "/model/shell/section_positions", {3U}, + section_positions.data(), section_positions.size(), + "ZETA", "1", "shell-natural", "section-position"); + { + Hdf5Handle dataset{ + RequireHdf5Id( + H5Dopen2(file, "/model/shell/section_positions", H5P_DEFAULT), + "Unable to reopen shell section positions."), + H5Dclose}; + WriteStringAttribute(dataset.Get(), "position_names", "BOTTOM,MIDDLE,TOP"); + } +} + +void WriteModel(const hid_t file, const Domain& domain, + const WriterModelData& model_data) { + WriteMetadata(file, domain); + WriteNodes(file, domain); + if (IsShellDomain(domain)) { + WriteShellElements(file, domain); + WriteShellModelData(file, domain, model_data); + } else { + WriteBeamElements(file, domain, model_data.beam_local_axes); + } +} + +} // namespace fesa::hdf5_internal diff --git a/src/fesa/io/hdf5/hdf5_model_writer.h b/src/fesa/io/hdf5/hdf5_model_writer.h new file mode 100644 index 0000000..18e57ee --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_model_writer.h @@ -0,0 +1,17 @@ +#ifndef FESA_IO_HDF5_HDF5_MODEL_WRITER_H_ +#define FESA_IO_HDF5_HDF5_MODEL_WRITER_H_ + +#include + +#include "fesa/model/domain.h" +#include "io/hdf5/hdf5_primitives.h" + +namespace fesa::hdf5_internal { + +/// @brief Writes metadata and exact model/source identity datasets. +void WriteModel(hid_t file, const Domain& domain, + const WriterModelData& model_data); + +} // namespace fesa::hdf5_internal + +#endif // FESA_IO_HDF5_HDF5_MODEL_WRITER_H_ diff --git a/src/fesa/io/hdf5/hdf5_primitives.cpp b/src/fesa/io/hdf5/hdf5_primitives.cpp new file mode 100644 index 0000000..c58076e --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_primitives.cpp @@ -0,0 +1,259 @@ +#include "io/hdf5/hdf5_primitives.h" + +#include +#include +#include +#include +#include +#include +#include + +namespace fesa::hdf5_internal { + +Status OutputFailure(const std::string& code, const std::string& message) { + return Status::Failure(FailureCategory::kOutput, + {{Severity::kError, code, {}, "", "", message}}); +} + +bool IsValidUtf8(const std::string& value) { + const auto* bytes = reinterpret_cast(value.data()); + std::size_t index = 0U; + while (index < value.size()) { + const unsigned char first = bytes[index]; + if (first <= 0x7fU) { + ++index; + continue; + } + + std::size_t continuation_count = 0U; + std::uint32_t code_point = 0U; + if (first >= 0xc2U && first <= 0xdfU) { + continuation_count = 1U; + code_point = first & 0x1fU; + } else if (first >= 0xe0U && first <= 0xefU) { + continuation_count = 2U; + code_point = first & 0x0fU; + } else if (first >= 0xf0U && first <= 0xf4U) { + continuation_count = 3U; + code_point = first & 0x07U; + } else { + return false; + } + if (index + continuation_count >= value.size()) { + return false; + } + for (std::size_t offset = 1U; offset <= continuation_count; ++offset) { + const unsigned char next = bytes[index + offset]; + if ((next & 0xc0U) != 0x80U) { + return false; + } + code_point = (code_point << 6U) | (next & 0x3fU); + } + if ((continuation_count == 2U && code_point < 0x800U) || + (continuation_count == 3U && code_point < 0x10000U) || + (code_point >= 0xd800U && code_point <= 0xdfffU) || + code_point > 0x10ffffU) { + return false; + } + index += continuation_count + 1U; + } + 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 IsShellDomain(const Domain& domain) noexcept { + return !domain.ShellElements().Empty(); +} + +const char* ShellSourceTypeName(const ShellSourceElementType type) { + return type == ShellSourceElementType::kS4 ? "S4" : "S4R"; +} + +bool SizeProductFits(const std::size_t left, const std::size_t right, + std::size_t& product) { + if (right != 0U && left > (std::numeric_limits::max)() / right) { + return false; + } + product = left * right; + return true; +} + +std::string NormalizedPathString(const std::filesystem::path& path) { + if (path.empty()) { + return {}; + } + return std::filesystem::absolute(path).lexically_normal().generic_u8string(); +} + +std::string SourceInputIdentity(const Domain& domain) { + return "path=" + NormalizedPathString(domain.SourcePath()) + + ";content_identity=" + domain.SourceContentIdentity(); +} + +Hdf5Handle MakeUtf8StringType() { + Hdf5Handle type{ + RequireHdf5Id(H5Tcopy(H5T_C_S1), "Unable to copy an HDF5 string type."), + H5Tclose}; + RequireHdf5(H5Tset_size(type.Get(), H5T_VARIABLE), + "Unable to configure a variable-length HDF5 string."); + RequireHdf5(H5Tset_cset(type.Get(), H5T_CSET_UTF8), + "Unable to configure UTF-8 HDF5 strings."); + RequireHdf5(H5Tset_strpad(type.Get(), H5T_STR_NULLTERM), + "Unable to configure HDF5 string padding."); + return type; +} + +void WriteStringAttribute(const hid_t object, const char* name, + const std::string& value) { + auto type = MakeUtf8StringType(); + Hdf5Handle space{RequireHdf5Id(H5Screate(H5S_SCALAR), + "Unable to create an attribute space."), + H5Sclose}; + Hdf5Handle attribute{ + RequireHdf5Id(H5Acreate2(object, name, type.Get(), space.Get(), + H5P_DEFAULT, H5P_DEFAULT), + "Unable to create an HDF5 string attribute."), + H5Aclose}; + const char* raw = value.c_str(); + RequireHdf5(H5Awrite(attribute.Get(), type.Get(), &raw), + "Unable to write an HDF5 string attribute."); +} + +void WriteUint64Attribute(const hid_t object, const char* name, + const std::uint64_t value) { + Hdf5Handle space{RequireHdf5Id(H5Screate(H5S_SCALAR), + "Unable to create an attribute space."), + H5Sclose}; + Hdf5Handle attribute{ + RequireHdf5Id(H5Acreate2(object, name, H5T_STD_U64LE, space.Get(), + H5P_DEFAULT, H5P_DEFAULT), + "Unable to create an HDF5 integer attribute."), + H5Aclose}; + RequireHdf5(H5Awrite(attribute.Get(), H5T_NATIVE_UINT64, &value), + "Unable to write an HDF5 integer attribute."); +} + +Hdf5Handle CreateGroup(const hid_t file, const char* path) { + Hdf5Handle link_properties{RequireHdf5Id(H5Pcreate(H5P_LINK_CREATE), + "Unable to create link properties."), + H5Pclose}; + RequireHdf5(H5Pset_create_intermediate_group(link_properties.Get(), 1U), + "Unable to enable intermediate HDF5 groups."); + return Hdf5Handle{RequireHdf5Id(H5Gcreate2(file, path, link_properties.Get(), + H5P_DEFAULT, H5P_DEFAULT), + "Unable to create an HDF5 group."), + H5Gclose}; +} + +Hdf5Handle CreateDatasetSpace(const std::vector& dimensions) { + return Hdf5Handle{ + RequireHdf5Id(H5Screate_simple(static_cast(dimensions.size()), + dimensions.data(), nullptr), + "Unable to create an HDF5 dataset space."), + H5Sclose}; +} + +void WriteResultAttributes(const hid_t dataset, + const std::string& component_names, + const std::string& component_units, + const std::string& coordinate_system, + const std::string& location) { + WriteStringAttribute(dataset, "component_names", component_names); + WriteStringAttribute(dataset, "component_unit_dimensions", component_units); + WriteStringAttribute(dataset, "coordinate_system", coordinate_system); + WriteStringAttribute(dataset, "location", location); + WriteStringAttribute(dataset, "step_name", kStepName); + WriteUint64Attribute(dataset, "frame_index", 0U); +} + +Hdf5Handle WriteDoubleValues(const hid_t file, const std::string& path, + const std::vector& dimensions, + const double* values, + const std::size_t value_count) { + auto space = CreateDatasetSpace(dimensions); + Hdf5Handle dataset{ + RequireHdf5Id(H5Dcreate2(file, path.c_str(), H5T_IEEE_F64LE, space.Get(), + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT), + "Unable to create a floating-point HDF5 dataset."), + H5Dclose}; + if (value_count != 0U) { + RequireHdf5(H5Dwrite(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values), + "Unable to write a floating-point HDF5 dataset."); + } + return dataset; +} + +void WriteDoubleDataset(const hid_t file, const std::string& path, + const std::vector& 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 = WriteDoubleValues(file, path, dimensions, values, value_count); + WriteResultAttributes(dataset.Get(), component_names, component_units, + coordinate_system, location); +} + +void WriteModelDoubleDataset(const hid_t file, const std::string& path, + const std::vector& 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 = WriteDoubleValues(file, path, dimensions, values, value_count); + WriteStringAttribute(dataset.Get(), "component_names", component_names); + WriteStringAttribute(dataset.Get(), "component_unit_dimensions", + component_units); + WriteStringAttribute(dataset.Get(), "coordinate_system", coordinate_system); + WriteStringAttribute(dataset.Get(), "location", location); +} + +void WriteUint8Dataset(const hid_t file, const std::string& path, + const std::vector& dimensions, + const std::uint8_t* values, + const std::size_t value_count) { + auto space = CreateDatasetSpace(dimensions); + Hdf5Handle dataset{ + RequireHdf5Id(H5Dcreate2(file, path.c_str(), H5T_STD_U8LE, space.Get(), + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT), + "Unable to create a uint8 HDF5 dataset."), + H5Dclose}; + if (value_count != 0U) { + RequireHdf5(H5Dwrite(dataset.Get(), H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values), + "Unable to write a uint8 HDF5 dataset."); + } + WriteStringAttribute(dataset.Get(), "component_names", + "UX,UY,UZ,URX,URY,URZ"); + WriteStringAttribute(dataset.Get(), "value_meaning", "0=free,1=constrained"); +} + +Hdf5Handle WriteCompoundDataset(const hid_t file, const char* path, + const std::size_t row_count, + const hid_t file_type, const hid_t memory_type, + const void* rows) { + const std::vector dimensions = {static_cast(row_count)}; + auto space = CreateDatasetSpace(dimensions); + Hdf5Handle dataset{ + RequireHdf5Id(H5Dcreate2(file, path, file_type, space.Get(), H5P_DEFAULT, + H5P_DEFAULT, H5P_DEFAULT), + "Unable to create a compound HDF5 dataset."), + H5Dclose}; + if (row_count != 0U) { + RequireHdf5(H5Dwrite(dataset.Get(), memory_type, H5S_ALL, H5S_ALL, + H5P_DEFAULT, rows), + "Unable to write a compound HDF5 dataset."); + } + return dataset; +} + +} // namespace fesa::hdf5_internal diff --git a/src/fesa/io/hdf5/hdf5_primitives.h b/src/fesa/io/hdf5/hdf5_primitives.h new file mode 100644 index 0000000..429448b --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_primitives.h @@ -0,0 +1,132 @@ +#ifndef FESA_IO_HDF5_HDF5_PRIMITIVES_H_ +#define FESA_IO_HDF5_HDF5_PRIMITIVES_H_ + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/core/source_identity.h" +#include "fesa/core/status.h" +#include "fesa/model/domain.h" +#include "io/hdf5/hdf5_raii.h" + +namespace fesa::hdf5_internal { + +inline constexpr std::size_t kDofsPerNode = 6U; +inline constexpr std::size_t kEndpointCount = 2U; +inline constexpr std::size_t kGaussPointCount = 2U; +inline constexpr std::size_t kEndActionComponentCount = 6U; +inline constexpr std::size_t kGeneralizedComponentCount = 4U; +inline constexpr std::size_t kShellNodeCount = 4U; +inline constexpr std::size_t kShellLocationCount = 4U; +inline constexpr std::size_t kShellGeneralizedComponentCount = 8U; +inline constexpr std::size_t kShellSectionPositionCount = 3U; +inline constexpr std::size_t kShellStressComponentCount = 3U; +inline constexpr const char* kStepName = "Step-1"; +inline constexpr std::size_t kFrameIndex = 0U; +inline constexpr const char* kStepRoot = "/steps/Step-1/frames/0"; + +using AxisSet = std::array; + +/// @brief Owns model values derived before candidate-file creation. +struct WriterModelData { + std::vector beam_local_axes; + std::vector constraint_mask; + std::vector prescribed_displacement; +}; + +/// @brief Describes one backend-neutral float64 result dataset. +struct DoubleDatasetPlan { + std::string path; + std::vector dimensions; + std::vector values; + std::string component_names; + std::string component_units; + std::string coordinate_system; + std::string location; +}; + +/// @brief Owns the fully validated model and result write plan. +struct WriterPlan { + WriterModelData model_data; + std::vector result_datasets; +}; + +/// @brief Creates one structured output failure. +Status OutputFailure(const std::string& code, const std::string& message); + +/// @brief Tests exact preserved source identity. +bool SameIdentity(const SourceEntityId& left, const SourceEntityId& right); + +/// @brief Tests a fixed component record for finite values. +template +bool IsFinite(const std::array& values) { + for (const double value : values) { + if (!std::isfinite(value)) { + return false; + } + } + return true; +} + +/// @brief Validates UTF-8 without accepting overlong or surrogate encodings. +bool IsValidUtf8(const std::string& value); + +/// @brief Identifies the schema-v0 shell variant. +bool IsShellDomain(const Domain& domain) noexcept; + +/// @brief Returns the preserved shell source element type spelling. +const char* ShellSourceTypeName(ShellSourceElementType type); + +/// @brief Computes a checked size product. +bool SizeProductFits(std::size_t left, std::size_t right, std::size_t& product); + +/// @brief Returns a normalized absolute UTF-8 path string. +std::string NormalizedPathString(const std::filesystem::path& path); + +/// @brief Returns schema-v0 source path and content identity metadata. +std::string SourceInputIdentity(const Domain& domain); + +Hdf5Handle MakeUtf8StringType(); +void WriteStringAttribute(hid_t object, const char* name, + const std::string& value); +void WriteUint64Attribute(hid_t object, const char* name, std::uint64_t value); +Hdf5Handle CreateGroup(hid_t file, const char* path); +Hdf5Handle CreateDatasetSpace(const std::vector& dimensions); +void WriteResultAttributes(hid_t dataset, const std::string& component_names, + const std::string& component_units, + const std::string& coordinate_system, + const std::string& location); +Hdf5Handle WriteDoubleValues(hid_t file, const std::string& path, + const std::vector& dimensions, + const double* values, std::size_t value_count); +void WriteDoubleDataset(hid_t file, const std::string& path, + const std::vector& dimensions, + const double* values, std::size_t value_count, + const std::string& component_names, + const std::string& component_units, + const std::string& coordinate_system, + const std::string& location); +void WriteModelDoubleDataset(hid_t file, const std::string& path, + const std::vector& dimensions, + const double* values, std::size_t value_count, + const std::string& component_names, + const std::string& component_units, + const std::string& coordinate_system, + const std::string& location); +void WriteUint8Dataset(hid_t file, const std::string& path, + const std::vector& dimensions, + const std::uint8_t* values, std::size_t value_count); +Hdf5Handle WriteCompoundDataset(hid_t file, const char* path, + std::size_t row_count, hid_t file_type, + hid_t memory_type, const void* rows); + +} // namespace fesa::hdf5_internal + +#endif // FESA_IO_HDF5_HDF5_PRIMITIVES_H_ diff --git a/src/fesa/io/hdf5/hdf5_raii.h b/src/fesa/io/hdf5/hdf5_raii.h new file mode 100644 index 0000000..051c473 --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_raii.h @@ -0,0 +1,110 @@ +#ifndef FESA_IO_HDF5_HDF5_RAII_H_ +#define FESA_IO_HDF5_HDF5_RAII_H_ + +#include + +#include +#include + +namespace fesa::hdf5_internal { + +/// @brief Owns one HDF5 identifier and its matching close operation. +/// @note This move-only wrapper never crosses the private HDF5 boundary. +class Hdf5Handle { + public: + using Closer = herr_t (*)(hid_t); + + Hdf5Handle() = default; + Hdf5Handle(hid_t value, Closer closer) : value_{value}, closer_{closer} {} + Hdf5Handle(const Hdf5Handle&) = delete; + Hdf5Handle& operator=(const Hdf5Handle&) = delete; + Hdf5Handle(Hdf5Handle&& other) noexcept + : value_{other.value_}, closer_{other.closer_} { + other.value_ = -1; + other.closer_ = nullptr; + } + Hdf5Handle& operator=(Hdf5Handle&& other) noexcept { + if (this != &other) { + Reset(); + value_ = other.value_; + closer_ = other.closer_; + other.value_ = -1; + other.closer_ = nullptr; + } + return *this; + } + ~Hdf5Handle() { Reset(); } + + hid_t Get() const noexcept { return value_; } + herr_t CloseChecked() noexcept { + if (value_ < 0 || closer_ == nullptr) { + return 0; + } + const hid_t value = value_; + const Closer closer = closer_; + value_ = -1; + closer_ = nullptr; + return closer(value); + } + + private: + void Reset() noexcept { + if (value_ >= 0 && closer_ != nullptr) { + (void)closer_(value_); + } + value_ = -1; + closer_ = nullptr; + } + + hid_t value_{-1}; + Closer closer_{nullptr}; +}; + +/// @brief Converts an HDF5 failure into the private exception transport. +class Hdf5Failure final : public std::runtime_error { + public: + explicit Hdf5Failure(const std::string& message) + : std::runtime_error{message} {} +}; + +/// @brief Requires a successful HDF5 status code. +inline void RequireHdf5(herr_t result, const char* message) { + if (result < 0) { + throw Hdf5Failure{message}; + } +} + +/// @brief Requires a valid HDF5 identifier. +inline hid_t RequireHdf5Id(hid_t result, const char* message) { + if (result < 0) { + throw Hdf5Failure{message}; + } + return result; +} + +/// @brief Suppresses backend error-stack output while FESA translates errors. +class Hdf5ErrorSilencer { + public: + Hdf5ErrorSilencer() { + if (H5Eget_auto2(H5E_DEFAULT, &callback_, &client_data_) >= 0 && + H5Eset_auto2(H5E_DEFAULT, nullptr, nullptr) >= 0) { + active_ = true; + } + } + Hdf5ErrorSilencer(const Hdf5ErrorSilencer&) = delete; + Hdf5ErrorSilencer& operator=(const Hdf5ErrorSilencer&) = delete; + ~Hdf5ErrorSilencer() { + if (active_) { + (void)H5Eset_auto2(H5E_DEFAULT, callback_, client_data_); + } + } + + private: + H5E_auto2_t callback_{nullptr}; + void* client_data_{nullptr}; + bool active_{false}; +}; + +} // namespace fesa::hdf5_internal + +#endif // FESA_IO_HDF5_HDF5_RAII_H_ diff --git a/src/fesa/io/hdf5/hdf5_result_writer.cpp b/src/fesa/io/hdf5/hdf5_result_writer.cpp new file mode 100644 index 0000000..779b839 --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_result_writer.cpp @@ -0,0 +1,955 @@ +#include "io/hdf5/hdf5_result_writer.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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, 3>& frame) { + constexpr double kTolerance = 1.0e-12; + const std::array 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, 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 locations{ + ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2, + ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4}; + const std::array, kShellLocationCount> coordinates{ + {{-gauss, -gauss}, {gauss, -gauss}, {gauss, gauss}, {-gauss, gauss}}}; + const std::array positions{ + ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle, + ShellSectionPosition::kTop}; + constexpr std::array 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& 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 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(row_index / kEndpointCount); + const int expected_endpoint = static_cast(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( + 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(row_index / kGaussPointCount); + const int expected_gauss_point = + static_cast(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(element_index) || + row.gauss_point != static_cast(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 FlattenEndpointValues( + const std::vector& rows, const bool section_resultants) { + const std::size_t components = section_resultants ? kGeneralizedComponentCount + : kEndActionComponentCount; + std::vector 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 FlattenGaussValues(const std::vector& rows, + const bool resultants) { + std::vector 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 rows; + rows.reserve(state.StressResults().size()); + for (const auto& row : state.StressResults()) { + rows.push_back({static_cast(row.element), + static_cast(row.gauss_point), + static_cast(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 local_frames; + std::vector generalized_strains; + std::vector section_resultants; + std::vector 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(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& input_diagnostics) { + std::vector diagnostics = input_diagnostics; + SortDiagnostics(diagnostics); + std::vector files; + files.reserve(diagnostics.size()); + for (const auto& diagnostic : diagnostics) { + files.push_back(NormalizedPathString(diagnostic.location.file)); + } + std::vector 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(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 nodal_dimensions = { + static_cast(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 endpoint_action_dimensions = { + static_cast(domain.BeamElements().Size()), kEndpointCount, + kEndActionComponentCount}; + const std::vector generalized_dimensions = { + static_cast(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 CopyVector(const Vector& source) { + std::vector values(source.Size()); + for (std::size_t index = 0U; index < source.Size(); ++index) { + values[index] = source[index]; + } + return values; +} + +std::vector BuildResultDatasetPlans( + const Domain& domain, const AnalysisState& state) { + std::vector datasets; + const std::vector 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 local_frames; + std::vector generalized_strains; + std::vector section_resultants; + std::vector 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(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(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& datasets) { + std::set 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 BuildWriterPlan(const std::filesystem::path& output_path, + const Domain& domain, + const AnalysisState& state, + const std::vector& diagnostics) { + WriterPlan plan; + const Status validation = ValidateWriterInput(output_path, domain, state, + diagnostics, plan.model_data); + if (!validation.IsOk()) { + return Result::Failure(validation); + } + plan.result_datasets = BuildResultDatasetPlans(domain, state); + const Status inventory = ValidateFiniteInventory(plan.result_datasets); + if (!inventory.IsOk()) { + return Result::Failure(inventory); + } + return Result::Success(std::move(plan)); +} + +void WriteResults(const hid_t file, const Domain& domain, + const AnalysisState& state, + const std::vector& 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 diff --git a/src/fesa/io/hdf5/hdf5_result_writer.h b/src/fesa/io/hdf5/hdf5_result_writer.h new file mode 100644 index 0000000..802dc70 --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_result_writer.h @@ -0,0 +1,33 @@ +#ifndef FESA_IO_HDF5_HDF5_RESULT_WRITER_H_ +#define FESA_IO_HDF5_HDF5_RESULT_WRITER_H_ + +#include + +#include +#include + +#include "fesa/analysis/analysis_state.h" +#include "fesa/core/diagnostic.h" +#include "fesa/core/status.h" +#include "fesa/model/domain.h" +#include "io/hdf5/hdf5_primitives.h" + +namespace fesa::hdf5_internal { + +/// @brief Builds the complete backend-neutral result dataset inventory. +Result BuildWriterPlan(const std::filesystem::path& output_path, + const Domain& domain, + const AnalysisState& state, + const std::vector& diagnostics); + +/// @brief Rejects malformed, duplicate, shape-mismatched, or nonfinite plans. +Status ValidateFiniteInventory(const std::vector& datasets); + +/// @brief Writes mandatory beam/shell/global results and diagnostics. +void WriteResults(hid_t file, const Domain& domain, const AnalysisState& state, + const std::vector& diagnostics, + const WriterPlan& plan); + +} // namespace fesa::hdf5_internal + +#endif // FESA_IO_HDF5_HDF5_RESULT_WRITER_H_ diff --git a/src/fesa/io/hdf5/hdf5_results_writer.cpp b/src/fesa/io/hdf5/hdf5_results_writer.cpp index f968da6..077ca81 100644 --- a/src/fesa/io/hdf5/hdf5_results_writer.cpp +++ b/src/fesa/io/hdf5/hdf5_results_writer.cpp @@ -1,2319 +1,73 @@ -#define NOMINMAX #include "fesa/io/hdf5/hdf5_results_writer.h" -#include -#include - -#include -#include -#include -#include -#include -#include +#include #include -#include -#include -#include #include #include -#include "fesa/analysis/analysis_model.h" -#include "fesa/build_info.h" -#include "fesa/fem/dof_manager.h" -#include "fesa/math/vector3.h" +#include "io/hdf5/hdf5_atomic_file.h" +#include "io/hdf5/hdf5_model_writer.h" +#include "io/hdf5/hdf5_primitives.h" +#include "io/hdf5/hdf5_raii.h" +#include "io/hdf5/hdf5_result_writer.h" +#include "io/hdf5/hdf5_self_check.h" namespace fesa { -namespace { - -constexpr std::size_t kDofsPerNode = 6U; -constexpr std::size_t kEndpointCount = 2U; -constexpr std::size_t kGaussPointCount = 2U; -constexpr std::size_t kEndActionComponentCount = 6U; -constexpr std::size_t kGeneralizedComponentCount = 4U; -constexpr std::size_t kShellNodeCount = 4U; -constexpr std::size_t kShellLocationCount = 4U; -constexpr std::size_t kShellGeneralizedComponentCount = 8U; -constexpr std::size_t kShellSectionPositionCount = 3U; -constexpr std::size_t kShellStressComponentCount = 3U; -constexpr const char* kStepName = "Step-1"; -constexpr std::size_t kFrameIndex = 0U; -constexpr const char* kStepRoot = "/steps/Step-1/frames/0"; - -/// @brief Owns one HDF5 backend identifier and its matching close operation. -/// @note Move-only lifetime keeps HDF5 handles inside the writer -/// implementation. -class Hdf5Handle { - public: - using Closer = herr_t (*)(hid_t); - - Hdf5Handle() = default; - Hdf5Handle(const hid_t value, Closer closer) - : value_{value}, closer_{closer} {} - Hdf5Handle(const Hdf5Handle&) = delete; - Hdf5Handle& operator=(const Hdf5Handle&) = delete; - Hdf5Handle(Hdf5Handle&& other) noexcept - : value_{other.value_}, closer_{other.closer_} { - other.value_ = -1; - other.closer_ = nullptr; - } - Hdf5Handle& operator=(Hdf5Handle&& other) noexcept { - if (this != &other) { - Reset(); - value_ = other.value_; - closer_ = other.closer_; - other.value_ = -1; - other.closer_ = nullptr; - } - return *this; - } - ~Hdf5Handle() { Reset(); } - - hid_t Get() const noexcept { return value_; } - herr_t CloseChecked() noexcept { - if (value_ < 0 || closer_ == nullptr) { - return 0; - } - const hid_t value = value_; - const Closer closer = closer_; - value_ = -1; - closer_ = nullptr; - return closer(value); - } - - private: - void Reset() noexcept { - if (value_ >= 0 && closer_ != nullptr) { - (void)closer_(value_); - } - value_ = -1; - closer_ = nullptr; - } - - hid_t value_{-1}; - Closer closer_{nullptr}; -}; - -/// @brief Temporarily suppresses HDF5's process-global automatic error printer. -/// @note Expected backend failures become one structured FESA diagnostic, not -/// an HDF5 error-stack dump mixed into deterministic CLI output. -class Hdf5ErrorSilencer { - public: - Hdf5ErrorSilencer() { - if (H5Eget_auto2(H5E_DEFAULT, &callback_, &client_data_) >= 0 && - H5Eset_auto2(H5E_DEFAULT, nullptr, nullptr) >= 0) { - active_ = true; - } - } - Hdf5ErrorSilencer(const Hdf5ErrorSilencer&) = delete; - Hdf5ErrorSilencer& operator=(const Hdf5ErrorSilencer&) = delete; - ~Hdf5ErrorSilencer() { - if (active_) { - (void)H5Eset_auto2(H5E_DEFAULT, callback_, client_data_); - } - } - - private: - H5E_auto2_t callback_{nullptr}; - void* client_data_{nullptr}; - bool active_{false}; -}; - -/// @brief Removes an uncommitted temporary result file on every failure path. -/// @note Release is called only after the atomic finalization operation -/// succeeds. -class TemporaryFileGuard { - public: - explicit TemporaryFileGuard(std::filesystem::path path) - : path_{std::move(path)} {} - TemporaryFileGuard(const TemporaryFileGuard&) = delete; - TemporaryFileGuard& operator=(const TemporaryFileGuard&) = delete; - ~TemporaryFileGuard() { - if (active_) { - std::error_code ignored; - (void)std::filesystem::remove(path_, ignored); - } - } - void Release() noexcept { active_ = false; } - - private: - std::filesystem::path path_; - bool active_{true}; -}; - -class Hdf5Failure final : public std::runtime_error { - public: - explicit Hdf5Failure(const std::string& message) - : std::runtime_error{message} {} -}; - -void RequireHdf5(const herr_t result, const char* message) { - if (result < 0) { - throw Hdf5Failure{message}; - } -} - -hid_t RequireHdf5Id(const hid_t result, const char* message) { - if (result < 0) { - throw Hdf5Failure{message}; - } - return result; -} - -Status OutputFailure(const std::string& code, const std::string& message) { - return Status::Failure(FailureCategory::kOutput, - {{Severity::kError, code, {}, "", "", message}}); -} - -template -bool IsFinite(const std::array& values) { - return std::all_of(values.begin(), values.end(), - [](const double value) { return std::isfinite(value); }); -} - -bool IsValidUtf8(const std::string& value) { - const auto* bytes = reinterpret_cast(value.data()); - std::size_t index = 0U; - while (index < value.size()) { - const unsigned char first = bytes[index]; - if (first <= 0x7fU) { - ++index; - continue; - } - - std::size_t continuation_count = 0U; - std::uint32_t code_point = 0U; - if (first >= 0xc2U && first <= 0xdfU) { - continuation_count = 1U; - code_point = first & 0x1fU; - } else if (first >= 0xe0U && first <= 0xefU) { - continuation_count = 2U; - code_point = first & 0x0fU; - } else if (first >= 0xf0U && first <= 0xf4U) { - continuation_count = 3U; - code_point = first & 0x07U; - } else { - return false; - } - if (index + continuation_count >= value.size()) { - return false; - } - for (std::size_t offset = 1U; offset <= continuation_count; ++offset) { - const unsigned char next = bytes[index + offset]; - if ((next & 0xc0U) != 0x80U) { - return false; - } - code_point = (code_point << 6U) | (next & 0x3fU); - } - if ((continuation_count == 2U && code_point < 0x800U) || - (continuation_count == 3U && code_point < 0x10000U) || - (code_point >= 0xd800U && code_point <= 0xdfffU) || - code_point > 0x10ffffU) { - return false; - } - index += continuation_count + 1U; - } - 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; -} - -using AxisSet = std::array; - -struct WriterModelData { - std::vector beam_local_axes; - std::vector constraint_mask; - std::vector prescribed_displacement; -}; - -bool IsShellDomain(const Domain& domain) noexcept { - return !domain.ShellElements().Empty(); -} - -const char* ShellSourceTypeName(const ShellSourceElementType type) { - return type == ShellSourceElementType::kS4 ? "S4" : "S4R"; -} - -bool IsOrthonormalRightHanded( - const std::array, 3>& frame) { - constexpr double kTolerance = 1.0e-12; - const std::array 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); -} - -bool SizeProductFits(const std::size_t left, const std::size_t right, - std::size_t& product) { - if (right != 0U && left > (std::numeric_limits::max)() / right) { - return false; - } - product = left * right; - return true; -} - -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, 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 locations{ - ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2, - ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4}; - const std::array, kShellLocationCount> coordinates{ - {{-gauss, -gauss}, {gauss, -gauss}, {gauss, gauss}, {-gauss, gauss}}}; - const std::array positions{ - ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle, - ShellSectionPosition::kTop}; - constexpr std::array 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& 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 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(row_index / kEndpointCount); - const int expected_endpoint = static_cast(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( - 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(row_index / kGaussPointCount); - const int expected_gauss_point = - static_cast(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(element_index) || - row.gauss_point != static_cast(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::string NormalizedPathString(const std::filesystem::path& path) { - if (path.empty()) { - return {}; - } - return std::filesystem::absolute(path).lexically_normal().generic_u8string(); -} - -std::string SourceInputIdentity(const Domain& domain) { - return "path=" + NormalizedPathString(domain.SourcePath()) + - ";content_identity=" + domain.SourceContentIdentity(); -} - -Hdf5Handle MakeUtf8StringType() { - Hdf5Handle type{ - RequireHdf5Id(H5Tcopy(H5T_C_S1), "Unable to copy an HDF5 string type."), - H5Tclose}; - RequireHdf5(H5Tset_size(type.Get(), H5T_VARIABLE), - "Unable to configure a variable-length HDF5 string."); - RequireHdf5(H5Tset_cset(type.Get(), H5T_CSET_UTF8), - "Unable to configure UTF-8 HDF5 strings."); - RequireHdf5(H5Tset_strpad(type.Get(), H5T_STR_NULLTERM), - "Unable to configure HDF5 string padding."); - return type; -} - -void WriteStringAttribute(const hid_t object, const char* name, - const std::string& value) { - auto type = MakeUtf8StringType(); - Hdf5Handle space{RequireHdf5Id(H5Screate(H5S_SCALAR), - "Unable to create an attribute space."), - H5Sclose}; - Hdf5Handle attribute{ - RequireHdf5Id(H5Acreate2(object, name, type.Get(), space.Get(), - H5P_DEFAULT, H5P_DEFAULT), - "Unable to create an HDF5 string attribute."), - H5Aclose}; - const char* raw = value.c_str(); - RequireHdf5(H5Awrite(attribute.Get(), type.Get(), &raw), - "Unable to write an HDF5 string attribute."); -} - -void WriteUint64Attribute(const hid_t object, const char* name, - const std::uint64_t value) { - Hdf5Handle space{RequireHdf5Id(H5Screate(H5S_SCALAR), - "Unable to create an attribute space."), - H5Sclose}; - Hdf5Handle attribute{ - RequireHdf5Id(H5Acreate2(object, name, H5T_STD_U64LE, space.Get(), - H5P_DEFAULT, H5P_DEFAULT), - "Unable to create an HDF5 integer attribute."), - H5Aclose}; - RequireHdf5(H5Awrite(attribute.Get(), H5T_NATIVE_UINT64, &value), - "Unable to write an HDF5 integer attribute."); -} - -Hdf5Handle CreateGroup(const hid_t file, const char* path) { - Hdf5Handle link_properties{RequireHdf5Id(H5Pcreate(H5P_LINK_CREATE), - "Unable to create link properties."), - H5Pclose}; - RequireHdf5(H5Pset_create_intermediate_group(link_properties.Get(), 1U), - "Unable to enable intermediate HDF5 groups."); - return Hdf5Handle{RequireHdf5Id(H5Gcreate2(file, path, link_properties.Get(), - H5P_DEFAULT, H5P_DEFAULT), - "Unable to create an HDF5 group."), - H5Gclose}; -} - -Hdf5Handle CreateDatasetSpace(const std::vector& dimensions) { - return Hdf5Handle{ - RequireHdf5Id(H5Screate_simple(static_cast(dimensions.size()), - dimensions.data(), nullptr), - "Unable to create an HDF5 dataset space."), - H5Sclose}; -} - -void WriteResultAttributes(const hid_t dataset, - const std::string& component_names, - const std::string& component_units, - const std::string& coordinate_system, - const std::string& location) { - WriteStringAttribute(dataset, "component_names", component_names); - WriteStringAttribute(dataset, "component_unit_dimensions", component_units); - WriteStringAttribute(dataset, "coordinate_system", coordinate_system); - WriteStringAttribute(dataset, "location", location); - WriteStringAttribute(dataset, "step_name", kStepName); - WriteUint64Attribute(dataset, "frame_index", 0U); -} - -Hdf5Handle WriteDoubleValues(const hid_t file, const std::string& path, - const std::vector& dimensions, - const double* values, - const std::size_t value_count) { - auto space = CreateDatasetSpace(dimensions); - Hdf5Handle dataset{ - RequireHdf5Id(H5Dcreate2(file, path.c_str(), H5T_IEEE_F64LE, space.Get(), - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT), - "Unable to create a floating-point HDF5 dataset."), - H5Dclose}; - if (value_count != 0U) { - RequireHdf5(H5Dwrite(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, - H5P_DEFAULT, values), - "Unable to write a floating-point HDF5 dataset."); - } - return dataset; -} - -void WriteDoubleDataset(const hid_t file, const std::string& path, - const std::vector& 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 = WriteDoubleValues(file, path, dimensions, values, value_count); - WriteResultAttributes(dataset.Get(), component_names, component_units, - coordinate_system, location); -} - -void WriteModelDoubleDataset(const hid_t file, const std::string& path, - const std::vector& 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 = WriteDoubleValues(file, path, dimensions, values, value_count); - WriteStringAttribute(dataset.Get(), "component_names", component_names); - WriteStringAttribute(dataset.Get(), "component_unit_dimensions", - component_units); - WriteStringAttribute(dataset.Get(), "coordinate_system", coordinate_system); - WriteStringAttribute(dataset.Get(), "location", location); -} - -void WriteUint8Dataset(const hid_t file, const std::string& path, - const std::vector& dimensions, - const std::uint8_t* values, - const std::size_t value_count) { - auto space = CreateDatasetSpace(dimensions); - Hdf5Handle dataset{ - RequireHdf5Id(H5Dcreate2(file, path.c_str(), H5T_STD_U8LE, space.Get(), - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT), - "Unable to create a uint8 HDF5 dataset."), - H5Dclose}; - if (value_count != 0U) { - RequireHdf5(H5Dwrite(dataset.Get(), H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL, - H5P_DEFAULT, values), - "Unable to write a uint8 HDF5 dataset."); - } - WriteStringAttribute(dataset.Get(), "component_names", - "UX,UY,UZ,URX,URY,URZ"); - WriteStringAttribute(dataset.Get(), "value_meaning", "0=free,1=constrained"); -} - -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; - std::uint64_t node_internal_ids[2]; - double local_axes[9]; -}; - -struct ShellElementWriteRow { - 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 ShellMaterialWriteRow { - std::uint64_t internal_material_id; - const char* name; - double youngs_modulus; - double poisson_ratio; -}; - -struct ShellSectionWriteRow { - 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 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; -}; - -Hdf5Handle WriteCompoundDataset(const hid_t file, const char* path, - const std::size_t row_count, - const hid_t file_type, const hid_t memory_type, - const void* rows) { - const std::vector dimensions = {static_cast(row_count)}; - auto space = CreateDatasetSpace(dimensions); - Hdf5Handle dataset{ - RequireHdf5Id(H5Dcreate2(file, path, file_type, space.Get(), H5P_DEFAULT, - H5P_DEFAULT, H5P_DEFAULT), - "Unable to create a compound HDF5 dataset."), - H5Dclose}; - if (row_count != 0U) { - RequireHdf5(H5Dwrite(dataset.Get(), memory_type, H5S_ALL, H5S_ALL, - H5P_DEFAULT, rows), - "Unable to write a compound HDF5 dataset."); - } - return dataset; -} - -void WriteMetadata(const hid_t file, const Domain& domain) { - auto metadata = CreateGroup(file, "/metadata"); - WriteUint64Attribute(metadata.Get(), "schema_version", 0U); - WriteStringAttribute(metadata.Get(), "solver_version", - std::string{SolverVersion()}); - WriteStringAttribute(metadata.Get(), "source_input_identity", - SourceInputIdentity(domain)); - WriteStringAttribute(metadata.Get(), "unit_system_label", - "user-consistent-unspecified"); - if (IsShellDomain(domain)) { - WriteStringAttribute(metadata.Get(), "feature_id", - "linear-static-mitc4-shell"); - WriteStringAttribute( - metadata.Get(), "coordinate_convention", - "global-cartesian; shell-local=(e1,e2,e3); positive-thickness=+zeta"); - WriteStringAttribute(metadata.Get(), "internal_formulation", "FESA-MITC4"); - WriteStringAttribute(metadata.Get(), "integration_rule", - "2x2x2-gauss; mitc4-edge-midpoint-shear"); - } else { - WriteStringAttribute(metadata.Get(), "feature_id", - "linear-static-3d-euler-beam"); - WriteStringAttribute(metadata.Get(), "coordinate_convention", - "global-cartesian; beam-local=(t,n1,t-cross-n1)"); - WriteStringAttribute(metadata.Get(), "element_formulation", - "B33-3D-Euler-Bernoulli"); - } - WriteStringAttribute(metadata.Get(), "step_name", kStepName); - WriteUint64Attribute(metadata.Get(), "frame_index", 0U); -} - -void WriteNodes(const hid_t file, const Domain& domain) { - std::vector 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(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 = MakeUtf8StringType(); - const hsize_t coordinate_dimensions[] = {3U}; - Hdf5Handle file_coordinates{ - RequireHdf5Id(H5Tarray_create2(H5T_IEEE_F64LE, 1, coordinate_dimensions), - "Unable to create the node coordinate file type."), - H5Tclose}; - Hdf5Handle memory_coordinates{ - RequireHdf5Id( - H5Tarray_create2(H5T_NATIVE_DOUBLE, 1, coordinate_dimensions), - "Unable to create the node coordinate memory type."), - H5Tclose}; - Hdf5Handle file_type{ - RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)), - "Unable to create the node file type."), - H5Tclose}; - Hdf5Handle memory_type{ - RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)), - "Unable to create the node memory type."), - H5Tclose}; - RequireHdf5(H5Tinsert(file_type.Get(), "internal_node_id", - HOFFSET(NodeWriteRow, internal_node_id), H5T_STD_U64LE), - "Unable to define the node internal ID field."); - RequireHdf5( - H5Tinsert(file_type.Get(), "instance_name", - HOFFSET(NodeWriteRow, instance_name), string_type.Get()), - "Unable to define the node instance field."); - RequireHdf5(H5Tinsert(file_type.Get(), "source_label", - HOFFSET(NodeWriteRow, source_label), string_type.Get()), - "Unable to define the node source-label field."); - RequireHdf5( - H5Tinsert(file_type.Get(), "coordinates", - HOFFSET(NodeWriteRow, coordinates), file_coordinates.Get()), - "Unable to define the node coordinate field."); - RequireHdf5( - H5Tinsert(memory_type.Get(), "internal_node_id", - HOFFSET(NodeWriteRow, internal_node_id), H5T_NATIVE_UINT64), - "Unable to define the node internal ID memory field."); - RequireHdf5( - H5Tinsert(memory_type.Get(), "instance_name", - HOFFSET(NodeWriteRow, instance_name), string_type.Get()), - "Unable to define the node instance memory field."); - RequireHdf5(H5Tinsert(memory_type.Get(), "source_label", - HOFFSET(NodeWriteRow, source_label), string_type.Get()), - "Unable to define the node source-label memory field."); - RequireHdf5( - H5Tinsert(memory_type.Get(), "coordinates", - HOFFSET(NodeWriteRow, coordinates), memory_coordinates.Get()), - "Unable to define the node coordinate memory field."); - - auto dataset = - WriteCompoundDataset(file, "/model/nodes", rows.size(), file_type.Get(), - memory_type.Get(), rows.data()); - WriteStringAttribute(dataset.Get(), "coordinate_system", "global-cartesian"); - WriteStringAttribute(dataset.Get(), "units_label", "length"); -} - -void WriteBeamElements(const hid_t file, const Domain& domain, - const std::vector& axes) { - std::vector rows; - rows.reserve(domain.BeamElements().Size()); - for (std::size_t index = 0U; index < domain.BeamElements().Size(); ++index) { - const auto& element = domain.BeamElements()[index]; - ElementWriteRow row{static_cast(index), - element.source_id.instance_name.c_str(), - element.source_id.source_label_text.c_str(), - {static_cast(element.node_indices[0U]), - static_cast(element.node_indices[1U])}, - {}}; - std::copy(axes[index].begin(), axes[index].end(), row.local_axes); - rows.push_back(row); - } - - auto string_type = MakeUtf8StringType(); - const hsize_t node_dimensions[] = {2U}; - const hsize_t axis_dimensions[] = {3U, 3U}; - Hdf5Handle file_nodes{ - RequireHdf5Id(H5Tarray_create2(H5T_STD_U64LE, 1, node_dimensions), - "Unable to create the element connectivity file type."), - H5Tclose}; - Hdf5Handle memory_nodes{ - RequireHdf5Id(H5Tarray_create2(H5T_NATIVE_UINT64, 1, node_dimensions), - "Unable to create the element connectivity memory type."), - H5Tclose}; - Hdf5Handle file_axes{ - RequireHdf5Id(H5Tarray_create2(H5T_IEEE_F64LE, 2, axis_dimensions), - "Unable to create the local-axis file type."), - H5Tclose}; - Hdf5Handle memory_axes{ - RequireHdf5Id(H5Tarray_create2(H5T_NATIVE_DOUBLE, 2, axis_dimensions), - "Unable to create the local-axis memory type."), - H5Tclose}; - Hdf5Handle file_type{ - RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)), - "Unable to create the element file type."), - H5Tclose}; - Hdf5Handle memory_type{ - RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)), - "Unable to create the element memory type."), - H5Tclose}; - const auto insert_fields = [&](const hid_t type, const hid_t integer_type, - const hid_t node_type, const hid_t axes_type) { - RequireHdf5( - H5Tinsert(type, "internal_element_id", - HOFFSET(ElementWriteRow, internal_element_id), integer_type), - "Unable to define the element internal ID field."); - RequireHdf5( - H5Tinsert(type, "instance_name", - HOFFSET(ElementWriteRow, instance_name), string_type.Get()), - "Unable to define the element instance field."); - RequireHdf5( - H5Tinsert(type, "source_label", HOFFSET(ElementWriteRow, source_label), - string_type.Get()), - "Unable to define the element source-label field."); - RequireHdf5( - H5Tinsert(type, "node_internal_ids", - HOFFSET(ElementWriteRow, node_internal_ids), node_type), - "Unable to define the element connectivity field."); - RequireHdf5(H5Tinsert(type, "local_axes", - HOFFSET(ElementWriteRow, local_axes), axes_type), - "Unable to define the element local-axis field."); - }; - insert_fields(file_type.Get(), H5T_STD_U64LE, file_nodes.Get(), - file_axes.Get()); - insert_fields(memory_type.Get(), H5T_NATIVE_UINT64, memory_nodes.Get(), - memory_axes.Get()); - auto dataset = - WriteCompoundDataset(file, "/model/elements", rows.size(), - file_type.Get(), memory_type.Get(), rows.data()); - WriteStringAttribute(dataset.Get(), "formulation", "B33-3D-Euler-Bernoulli"); -} - -void WriteShellElements(const hid_t file, const Domain& domain) { - std::vector 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(index), - element.source_id.instance_name.c_str(), - element.source_id.source_label_text.c_str(), - ShellSourceTypeName(element.source_type), - kMitc4InternalFormulation.data(), - {static_cast(element.node_indices[0U]), - static_cast(element.node_indices[1U]), - static_cast(element.node_indices[2U]), - static_cast(element.node_indices[3U])}, - static_cast(element.section_index), - static_cast(element.material_index)}); - } - - auto string_type = MakeUtf8StringType(); - const hsize_t node_dimensions[] = {kShellNodeCount}; - Hdf5Handle file_nodes{ - RequireHdf5Id(H5Tarray_create2(H5T_STD_U64LE, 1, node_dimensions), - "Unable to create shell connectivity file type."), - H5Tclose}; - Hdf5Handle memory_nodes{ - RequireHdf5Id(H5Tarray_create2(H5T_NATIVE_UINT64, 1, node_dimensions), - "Unable to create shell connectivity memory type."), - H5Tclose}; - Hdf5Handle file_type{ - RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellElementWriteRow)), - "Unable to create shell element file type."), - H5Tclose}; - Hdf5Handle memory_type{ - RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellElementWriteRow)), - "Unable to create shell element memory type."), - H5Tclose}; - const auto insert_fields = [&](const hid_t type, const hid_t integer_type, - const hid_t node_type) { - RequireHdf5(H5Tinsert(type, "internal_element_id", - HOFFSET(ShellElementWriteRow, internal_element_id), - integer_type), - "Unable to define shell element ID field."); - RequireHdf5(H5Tinsert(type, "instance_name", - HOFFSET(ShellElementWriteRow, instance_name), - string_type.Get()), - "Unable to define shell element instance field."); - RequireHdf5(H5Tinsert(type, "source_label", - HOFFSET(ShellElementWriteRow, source_label), - string_type.Get()), - "Unable to define shell element source-label field."); - RequireHdf5(H5Tinsert(type, "source_element_type", - HOFFSET(ShellElementWriteRow, source_element_type), - string_type.Get()), - "Unable to define shell source-type field."); - RequireHdf5(H5Tinsert(type, "internal_formulation", - HOFFSET(ShellElementWriteRow, internal_formulation), - string_type.Get()), - "Unable to define shell formulation field."); - RequireHdf5( - H5Tinsert(type, "node_internal_ids", - HOFFSET(ShellElementWriteRow, node_internal_ids), node_type), - "Unable to define shell connectivity field."); - RequireHdf5( - H5Tinsert(type, "shell_section_internal_id", - HOFFSET(ShellElementWriteRow, shell_section_internal_id), - integer_type), - "Unable to define shell section ID field."); - RequireHdf5(H5Tinsert(type, "material_internal_id", - HOFFSET(ShellElementWriteRow, material_internal_id), - integer_type), - "Unable to define shell material ID field."); - }; - insert_fields(file_type.Get(), H5T_STD_U64LE, file_nodes.Get()); - insert_fields(memory_type.Get(), H5T_NATIVE_UINT64, memory_nodes.Get()); - auto dataset = - WriteCompoundDataset(file, "/model/elements", rows.size(), - file_type.Get(), memory_type.Get(), rows.data()); - WriteStringAttribute(dataset.Get(), "formulation", "FESA-MITC4"); -} - -void WriteShellMaterials(const hid_t file, const Domain& domain) { - std::vector 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(index), material.name.c_str(), - material.youngs_modulus, material.poisson_ratio}); - } - auto string_type = MakeUtf8StringType(); - Hdf5Handle file_type{ - RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellMaterialWriteRow)), - "Unable to create shell material file type."), - H5Tclose}; - Hdf5Handle memory_type{ - RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellMaterialWriteRow)), - "Unable to create shell material 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_material_id", - HOFFSET(ShellMaterialWriteRow, internal_material_id), - integer_type), - "Unable to define shell material ID field."); - RequireHdf5(H5Tinsert(type, "name", HOFFSET(ShellMaterialWriteRow, name), - string_type.Get()), - "Unable to define shell material name field."); - RequireHdf5( - H5Tinsert(type, "E", HOFFSET(ShellMaterialWriteRow, youngs_modulus), - double_type), - "Unable to define shell material E field."); - RequireHdf5( - H5Tinsert(type, "nu", HOFFSET(ShellMaterialWriteRow, poisson_ratio), - double_type), - "Unable to define shell material nu 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, "/model/shell/materials", rows.size(), - file_type.Get(), memory_type.Get(), rows.data()); - WriteStringAttribute(dataset.Get(), "component_unit_dimensions", - "force/length^2,1"); -} - -void WriteShellSections(const hid_t file, const Domain& domain) { - std::vector source_files; - source_files.reserve(domain.ShellSections().Size()); - for (std::size_t index = 0U; index < domain.ShellSections().Size(); ++index) { - const auto& section = domain.ShellSections()[index]; - source_files.push_back(NormalizedPathString(section.location.file)); - } - std::vector 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(index), - source_files[index].c_str(), - static_cast(section.location.line), - section.name.c_str(), - static_cast(section.material_index), - section.thickness}); - } - auto string_type = MakeUtf8StringType(); - Hdf5Handle file_type{ - RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellSectionWriteRow)), - "Unable to create shell section file type."), - H5Tclose}; - Hdf5Handle memory_type{ - RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(ShellSectionWriteRow)), - "Unable to create shell section 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_section_id", - HOFFSET(ShellSectionWriteRow, internal_section_id), - integer_type), - "Unable to define shell section ID field."); - RequireHdf5(H5Tinsert(type, "source_file", - HOFFSET(ShellSectionWriteRow, source_file), - string_type.Get()), - "Unable to define shell section source-file field."); - RequireHdf5( - H5Tinsert(type, "source_line", - HOFFSET(ShellSectionWriteRow, source_line), integer_type), - "Unable to define shell section source-line field."); - RequireHdf5(H5Tinsert(type, "source_elset", - HOFFSET(ShellSectionWriteRow, source_elset), - string_type.Get()), - "Unable to define shell section ELSET field."); - RequireHdf5(H5Tinsert(type, "material_internal_id", - HOFFSET(ShellSectionWriteRow, material_internal_id), - integer_type), - "Unable to define shell section material ID field."); - RequireHdf5( - H5Tinsert(type, "thickness", HOFFSET(ShellSectionWriteRow, thickness), - double_type), - "Unable to define shell section thickness 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, "/model/shell/sections", rows.size(), - file_type.Get(), memory_type.Get(), rows.data()); - WriteStringAttribute(dataset.Get(), "thickness_unit_dimension", "length"); - WriteStringAttribute(dataset.Get(), "layering", "centered-single-layer"); -} - -void WriteShellModelData(const hid_t file, const Domain& domain, - const WriterModelData& model_data) { - std::vector directors; - directors.reserve(domain.Nodes().size() * 3U); - std::vector frames; - 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()); - } - WriteModelDoubleDataset(file, "/model/shell/nodal_director", - {static_cast(domain.Nodes().size()), 3U}, - directors.data(), directors.size(), "D1,D2,D3", - "1,1,1", "global-cartesian", "nodal"); - WriteModelDoubleDataset(file, "/model/shell/nodal_frame", - {static_cast(domain.Nodes().size()), 3U, 3U}, - frames.data(), frames.size(), "X,Y,Z", "1,1,1", - "global-cartesian", "nodal-frame"); - { - Hdf5Handle dataset{ - RequireHdf5Id(H5Dopen2(file, "/model/shell/nodal_frame", H5P_DEFAULT), - "Unable to reopen the nodal-frame dataset."), - H5Dclose}; - WriteStringAttribute(dataset.Get(), "axis_names", "A,B,D"); - } - - WriteShellMaterials(file, domain); - WriteShellSections(file, domain); - const std::vector nodal_dimensions{ - static_cast(domain.Nodes().size()), kDofsPerNode}; - WriteUint8Dataset(file, "/model/nodal_constraint_mask", nodal_dimensions, - model_data.constraint_mask.data(), - model_data.constraint_mask.size()); - WriteModelDoubleDataset( - file, "/model/prescribed_displacement", nodal_dimensions, - model_data.prescribed_displacement.data(), - model_data.prescribed_displacement.size(), "UX,UY,UZ,URX,URY,URZ", - "length,length,length,radian,radian,radian", "global-cartesian", - "nodal-prescribed-value"); - - const double gauss = 1.0 / std::sqrt(3.0); - const std::array locations{-gauss, -gauss, gauss, -gauss, - gauss, gauss, -gauss, gauss}; - WriteModelDoubleDataset(file, "/model/shell/midsurface_locations", {4U, 2U}, - locations.data(), locations.size(), "XI,ETA", "1,1", - "shell-natural", "midsurface-location"); - const std::array section_positions{-1.0, 0.0, 1.0}; - WriteModelDoubleDataset(file, "/model/shell/section_positions", {3U}, - section_positions.data(), section_positions.size(), - "ZETA", "1", "shell-natural", "section-position"); - { - Hdf5Handle dataset{ - RequireHdf5Id( - H5Dopen2(file, "/model/shell/section_positions", H5P_DEFAULT), - "Unable to reopen shell section positions."), - H5Dclose}; - WriteStringAttribute(dataset.Get(), "position_names", "BOTTOM,MIDDLE,TOP"); - } -} - -std::vector FlattenEndpointValues( - const std::vector& rows, const bool section_resultants) { - const std::size_t components = section_resultants ? kGeneralizedComponentCount - : kEndActionComponentCount; - std::vector 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 FlattenGaussValues(const std::vector& rows, - const bool resultants) { - std::vector 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 rows; - rows.reserve(state.StressResults().size()); - for (const auto& row : state.StressResults()) { - rows.push_back({static_cast(row.element), - static_cast(row.gauss_point), - static_cast(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 local_frames; - std::vector generalized_strains; - std::vector section_resultants; - std::vector 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(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& input_diagnostics) { - std::vector diagnostics = input_diagnostics; - SortDiagnostics(diagnostics); - std::vector files; - files.reserve(diagnostics.size()); - for (const auto& diagnostic : diagnostics) { - files.push_back(NormalizedPathString(diagnostic.location.file)); - } - std::vector 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(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 nodal_dimensions = { - static_cast(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 endpoint_action_dimensions = { - static_cast(domain.BeamElements().Size()), kEndpointCount, - kEndActionComponentCount}; - const std::vector generalized_dimensions = { - static_cast(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); -} - -/// @brief Writes the complete candidate schema to a same-directory temporary. -/// @note The HDF5 file handle is closed before self-check and finalization. -void WriteFile(const std::filesystem::path& temporary_path, - const Domain& domain, const AnalysisState& state, - const std::vector& diagnostics, - const WriterModelData& model_data) { - Hdf5Handle file{ - RequireHdf5Id(H5Fcreate(temporary_path.string().c_str(), H5F_ACC_EXCL, - H5P_DEFAULT, H5P_DEFAULT), - "Unable to create the temporary HDF5 file."), - H5Fclose}; - WriteMetadata(file.Get(), domain); - (void)CreateGroup(file.Get(), "/model"); - (void)CreateGroup(file.Get(), "/steps/Step-1/frames/0/nodal"); - (void)CreateGroup(file.Get(), "/steps/Step-1/frames/0/element"); - if (IsShellDomain(domain)) { - (void)CreateGroup(file.Get(), "/model/shell"); - (void)CreateGroup(file.Get(), "/steps/Step-1/frames/0/element/shell"); - (void)CreateGroup(file.Get(), "/steps/Step-1/frames/0/global"); - } - WriteNodes(file.Get(), domain); - if (IsShellDomain(domain)) { - WriteShellElements(file.Get(), domain); - WriteShellModelData(file.Get(), domain, model_data); - } else { - WriteBeamElements(file.Get(), domain, model_data.beam_local_axes); - } - WriteResultDatasets(file.Get(), domain, state); - WriteDiagnostics(file.Get(), diagnostics); - RequireHdf5(H5Fflush(file.Get(), H5F_SCOPE_GLOBAL), - "Unable to flush the temporary HDF5 file."); - RequireHdf5(file.CloseChecked(), - "Unable to close the temporary HDF5 file after writing."); -} - -Hdf5Handle OpenDatasetForCheck(const hid_t file, const char* path) { - return Hdf5Handle{RequireHdf5Id(H5Dopen2(file, path, H5P_DEFAULT), - "A mandatory HDF5 dataset is missing."), - H5Dclose}; -} - -std::vector CheckedDimensions(const hid_t dataset) { - Hdf5Handle space{RequireHdf5Id(H5Dget_space(dataset), - "Unable to inspect an HDF5 dataspace."), - H5Sclose}; - const int rank = H5Sget_simple_extent_ndims(space.Get()); - if (rank < 0) { - throw Hdf5Failure{"Unable to inspect an HDF5 dataset rank."}; - } - std::vector dimensions(static_cast(rank)); - if (rank > 0) { - RequireHdf5( - H5Sget_simple_extent_dims(space.Get(), dimensions.data(), nullptr), - "Unable to inspect HDF5 dataset dimensions."); - } - return dimensions; -} - -void RequireStringAttribute(const hid_t object, const char* name, - const std::string& expected) { - Hdf5Handle attribute{ - RequireHdf5Id(H5Aopen(object, name, H5P_DEFAULT), - "A mandatory HDF5 string attribute is missing."), - H5Aclose}; - Hdf5Handle type{RequireHdf5Id(H5Aget_type(attribute.Get()), - "Unable to inspect an HDF5 attribute type."), - H5Tclose}; - if (H5Tget_class(type.Get()) != H5T_STRING || - H5Tis_variable_str(type.Get()) <= 0 || - H5Tget_cset(type.Get()) != H5T_CSET_UTF8) { - throw Hdf5Failure{"An HDF5 string attribute is not variable-length UTF-8."}; - } - char* raw = nullptr; - RequireHdf5(H5Aread(attribute.Get(), type.Get(), &raw), - "Unable to read an HDF5 string attribute."); - const std::string actual = raw == nullptr ? std::string{} : std::string{raw}; - if (raw != nullptr) { - RequireHdf5(H5free_memory(raw), "Unable to release an HDF5 string value."); - } - if (actual != expected) { - throw Hdf5Failure{"An HDF5 string attribute has the wrong value."}; - } -} - -void RequireUint64Attribute(const hid_t object, const char* name, - const std::uint64_t expected) { - Hdf5Handle attribute{ - RequireHdf5Id(H5Aopen(object, name, H5P_DEFAULT), - "A mandatory HDF5 integer attribute is missing."), - H5Aclose}; - Hdf5Handle type{ - RequireHdf5Id(H5Aget_type(attribute.Get()), - "Unable to inspect an HDF5 integer attribute type."), - H5Tclose}; - if (H5Tget_class(type.Get()) != H5T_INTEGER || - H5Tget_size(type.Get()) != sizeof(std::uint64_t) || - H5Tget_sign(type.Get()) != H5T_SGN_NONE || - H5Tequal(type.Get(), H5T_STD_U64LE) <= 0) { - throw Hdf5Failure{"An HDF5 integer attribute is not portable uint64."}; - } - std::uint64_t actual = 0U; - RequireHdf5(H5Aread(attribute.Get(), H5T_NATIVE_UINT64, &actual), - "Unable to read an HDF5 integer attribute."); - if (actual != expected) { - throw Hdf5Failure{"An HDF5 integer attribute has the wrong value."}; - } -} - -void RequirePortableCompoundMember(const hid_t compound_type, - const unsigned index, - const std::string& name) { - Hdf5Handle member_type{ - RequireHdf5Id(H5Tget_member_type(compound_type, index), - "Unable to inspect a compound HDF5 field type."), - H5Tclose}; - const bool is_uint64 = - name == "internal_node_id" || name == "internal_element_id" || - name == "internal_material_id" || name == "internal_section_id" || - name == "shell_section_internal_id" || name == "material_internal_id" || - name == "source_line" || name == "gauss_point_index" || - name == "section_point_index" || name == "line"; - const bool is_float64 = name == "x1" || name == "x2" || name == "S11" || - name == "E" || name == "nu" || name == "thickness"; - const bool is_string = - name == "instance_name" || name == "source_label" || - name == "source_element_type" || name == "internal_formulation" || - name == "name" || name == "source_file" || name == "source_elset" || - name == "source" || name == "severity" || name == "code" || - name == "file" || name == "keyword" || name == "entity_identity" || - name == "message"; - if (is_uint64) { - if (H5Tget_class(member_type.Get()) != H5T_INTEGER || - H5Tget_size(member_type.Get()) != sizeof(std::uint64_t) || - H5Tget_sign(member_type.Get()) != H5T_SGN_NONE || - H5Tequal(member_type.Get(), H5T_STD_U64LE) <= 0) { - throw Hdf5Failure{"A compound HDF5 field is not portable uint64."}; - } - return; - } - if (is_float64) { - if (H5Tget_class(member_type.Get()) != H5T_FLOAT || - H5Tget_size(member_type.Get()) != sizeof(double) || - H5Tequal(member_type.Get(), H5T_IEEE_F64LE) <= 0) { - throw Hdf5Failure{"A compound HDF5 field is not float64."}; - } - return; - } - if (is_string) { - if (H5Tget_class(member_type.Get()) != H5T_STRING || - H5Tis_variable_str(member_type.Get()) <= 0 || - H5Tget_cset(member_type.Get()) != H5T_CSET_UTF8) { - throw Hdf5Failure{"A compound HDF5 field is not UTF-8."}; - } - return; - } - - if (H5Tget_class(member_type.Get()) != H5T_ARRAY) { - throw Hdf5Failure{"A compound HDF5 array field has the wrong type."}; - } - const int rank = H5Tget_array_ndims(member_type.Get()); - if (rank <= 0) { - throw Hdf5Failure{"Unable to inspect a compound HDF5 array rank."}; - } - std::vector dimensions(static_cast(rank)); - RequireHdf5(H5Tget_array_dims2(member_type.Get(), dimensions.data()), - "Unable to inspect compound HDF5 array dimensions."); - Hdf5Handle base_type{ - RequireHdf5Id(H5Tget_super(member_type.Get()), - "Unable to inspect a compound HDF5 array base type."), - H5Tclose}; - if (name == "node_internal_ids") { - if ((dimensions != std::vector{2U} && - dimensions != std::vector{4U}) || - H5Tequal(base_type.Get(), H5T_STD_U64LE) <= 0) { - throw Hdf5Failure{"Element connectivity is not uint64[2] or uint64[4]."}; - } - } else if (name == "coordinates") { - if (dimensions != std::vector{3U} || - H5Tequal(base_type.Get(), H5T_IEEE_F64LE) <= 0) { - throw Hdf5Failure{"Node coordinates are not float64[3]."}; - } - } else if (name == "local_axes") { - if (dimensions != std::vector{3U, 3U} || - H5Tequal(base_type.Get(), H5T_IEEE_F64LE) <= 0) { - throw Hdf5Failure{"Element local axes are not float64[3,3]."}; - } - } else { - throw Hdf5Failure{"A compound HDF5 field has an unknown type contract."}; - } -} - -void RequireResultAttributes(const hid_t dataset, const char* component_names, - const char* component_units, - const char* coordinate_system, - const char* location) { - RequireStringAttribute(dataset, "component_names", component_names); - RequireStringAttribute(dataset, "component_unit_dimensions", component_units); - RequireStringAttribute(dataset, "coordinate_system", coordinate_system); - RequireStringAttribute(dataset, "location", location); - RequireStringAttribute(dataset, "step_name", kStepName); - RequireUint64Attribute(dataset, "frame_index", 0U); -} - -void RequireShellResultIdentity(const hid_t file, const char* path, - const bool uses_section_positions) { - auto dataset = OpenDatasetForCheck(file, path); - RequireStringAttribute(dataset.Get(), "source_element_type_dataset", - "/model/elements.source_element_type"); - RequireStringAttribute(dataset.Get(), "internal_formulation", "FESA-MITC4"); - RequireStringAttribute(dataset.Get(), "midsurface_location_dataset", - "/model/shell/midsurface_locations"); - RequireStringAttribute(dataset.Get(), "local_frame_dataset", - "/steps/Step-1/frames/0/element/shell/local_frame"); - if (uses_section_positions) { - RequireStringAttribute(dataset.Get(), "section_position_dataset", - "/model/shell/section_positions"); - } -} - -void RequireDoubleDataset(const hid_t file, const char* path, - const std::vector& expected_dimensions, - const char* component_names, - const char* component_units, - const char* coordinate_system, const char* location) { - auto dataset = OpenDatasetForCheck(file, path); - if (CheckedDimensions(dataset.Get()) != expected_dimensions) { - throw Hdf5Failure{"A floating-point HDF5 dataset has the wrong shape."}; - } - Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), - "Unable to inspect an HDF5 dataset type."), - H5Tclose}; - if (H5Tget_class(type.Get()) != H5T_FLOAT || - H5Tget_size(type.Get()) != sizeof(double) || - H5Tequal(type.Get(), H5T_IEEE_F64LE) <= 0) { - throw Hdf5Failure{ - "A result dataset is not IEEE-754 float64 little-endian."}; - } - RequireResultAttributes(dataset.Get(), component_names, component_units, - coordinate_system, location); - std::size_t value_count = 1U; - for (const hsize_t dimension : expected_dimensions) { - std::size_t next = 0U; - if (!SizeProductFits(value_count, static_cast(dimension), - next)) { - throw Hdf5Failure{"A result dataset shape overflows size_t."}; - } - value_count = next; - } - std::vector values(value_count); - if (!values.empty()) { - RequireHdf5(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, - H5P_DEFAULT, values.data()), - "Unable to read a result dataset during self-check."); - } - if (!std::all_of(values.begin(), values.end(), - [](const double value) { return std::isfinite(value); })) { - throw Hdf5Failure{"A result dataset contains a nonfinite value."}; - } -} - -void RequireModelDoubleDataset( - const hid_t file, const char* path, - const std::vector& expected_dimensions, - const char* component_names, const char* component_units, - const char* coordinate_system, const char* location, - const std::vector* expected_values = nullptr) { - auto dataset = OpenDatasetForCheck(file, path); - if (CheckedDimensions(dataset.Get()) != expected_dimensions) { - throw Hdf5Failure{"A model floating-point dataset has the wrong shape."}; - } - Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), - "Unable to inspect a model dataset type."), - H5Tclose}; - if (H5Tget_class(type.Get()) != H5T_FLOAT || - H5Tget_size(type.Get()) != sizeof(double) || - H5Tequal(type.Get(), H5T_IEEE_F64LE) <= 0) { - throw Hdf5Failure{"A model dataset is not IEEE-754 float64 little-endian."}; - } - RequireStringAttribute(dataset.Get(), "component_names", component_names); - RequireStringAttribute(dataset.Get(), "component_unit_dimensions", - component_units); - RequireStringAttribute(dataset.Get(), "coordinate_system", coordinate_system); - RequireStringAttribute(dataset.Get(), "location", location); - std::size_t value_count = 1U; - for (const hsize_t dimension : expected_dimensions) { - std::size_t next = 0U; - if (!SizeProductFits(value_count, static_cast(dimension), - next)) { - throw Hdf5Failure{"A model dataset shape overflows size_t."}; - } - value_count = next; - } - std::vector values(value_count); - if (!values.empty()) { - RequireHdf5(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, - H5P_DEFAULT, values.data()), - "Unable to read a model dataset during self-check."); - } - if (!std::all_of(values.begin(), values.end(), - [](const double value) { return std::isfinite(value); })) { - throw Hdf5Failure{"A model dataset contains a nonfinite value."}; - } - if (expected_values != nullptr && values != *expected_values) { - throw Hdf5Failure{"A fixed model dataset has the wrong value order."}; - } -} - -void RequireUint8Dataset(const hid_t file, const char* path, - const std::vector& expected_dimensions, - const std::vector& expected_values) { - auto dataset = OpenDatasetForCheck(file, path); - if (CheckedDimensions(dataset.Get()) != expected_dimensions) { - throw Hdf5Failure{"A uint8 HDF5 dataset has the wrong shape."}; - } - Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), - "Unable to inspect a uint8 dataset type."), - H5Tclose}; - if (H5Tget_class(type.Get()) != H5T_INTEGER || - H5Tget_size(type.Get()) != sizeof(std::uint8_t) || - H5Tget_sign(type.Get()) != H5T_SGN_NONE || - H5Tequal(type.Get(), H5T_STD_U8LE) <= 0) { - throw Hdf5Failure{"A constraint mask is not portable uint8."}; - } - std::vector values(expected_values.size()); - if (!values.empty()) { - RequireHdf5(H5Dread(dataset.Get(), H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL, - H5P_DEFAULT, values.data()), - "Unable to read the constraint mask during self-check."); - } - if (values != expected_values) { - throw Hdf5Failure{"The constraint mask has the wrong value order."}; - } - RequireStringAttribute(dataset.Get(), "component_names", - "UX,UY,UZ,URX,URY,URZ"); - RequireStringAttribute(dataset.Get(), "value_meaning", - "0=free,1=constrained"); -} - -void RequireCompoundDataset(const hid_t file, const char* path, - const hsize_t expected_rows, - const std::vector& expected_members) { - auto dataset = OpenDatasetForCheck(file, path); - if (CheckedDimensions(dataset.Get()) != std::vector{expected_rows}) { - throw Hdf5Failure{"A compound HDF5 dataset has the wrong shape."}; - } - Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), - "Unable to inspect a compound type."), - H5Tclose}; - if (H5Tget_class(type.Get()) != H5T_COMPOUND || - H5Tget_nmembers(type.Get()) != - static_cast(expected_members.size())) { - throw Hdf5Failure{"A compound HDF5 dataset has the wrong field count."}; - } - for (std::size_t index = 0U; index < expected_members.size(); ++index) { - char* raw_name = - H5Tget_member_name(type.Get(), static_cast(index)); - if (raw_name == nullptr) { - throw Hdf5Failure{"Unable to inspect a compound HDF5 field name."}; - } - const std::string actual_name{raw_name}; - RequireHdf5(H5free_memory(raw_name), - "Unable to release a compound field name."); - if (actual_name != expected_members[index]) { - throw Hdf5Failure{"A compound HDF5 field has the wrong identity."}; - } - RequirePortableCompoundMember(type.Get(), static_cast(index), - expected_members[index]); - } -} - -/// @brief Reopens the closed candidate and verifies its required external -/// schema. -void SelfCheckFile(const std::filesystem::path& path, const Domain& domain, - const AnalysisState& state, - const std::size_t diagnostic_count, - const WriterModelData& model_data) { - if (H5Fis_hdf5(path.string().c_str()) <= 0) { - throw Hdf5Failure{"The temporary output is not an HDF5 file."}; - } - Hdf5Handle file{ - RequireHdf5Id(H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), - "Unable to reopen the temporary HDF5 file read-only."), - H5Fclose}; - Hdf5Handle metadata{ - RequireHdf5Id(H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), - "The HDF5 metadata group is missing."), - H5Gclose}; - RequireUint64Attribute(metadata.Get(), "schema_version", 0U); - RequireStringAttribute(metadata.Get(), "solver_version", - std::string{SolverVersion()}); - RequireStringAttribute(metadata.Get(), "source_input_identity", - SourceInputIdentity(domain)); - RequireStringAttribute(metadata.Get(), "unit_system_label", - "user-consistent-unspecified"); - if (IsShellDomain(domain)) { - RequireStringAttribute(metadata.Get(), "feature_id", - "linear-static-mitc4-shell"); - RequireStringAttribute( - metadata.Get(), "coordinate_convention", - "global-cartesian; shell-local=(e1,e2,e3); positive-thickness=+zeta"); - RequireStringAttribute(metadata.Get(), "internal_formulation", - "FESA-MITC4"); - RequireStringAttribute(metadata.Get(), "integration_rule", - "2x2x2-gauss; mitc4-edge-midpoint-shear"); - } else { - RequireStringAttribute(metadata.Get(), "feature_id", - "linear-static-3d-euler-beam"); - RequireStringAttribute(metadata.Get(), "coordinate_convention", - "global-cartesian; beam-local=(t,n1,t-cross-n1)"); - RequireStringAttribute(metadata.Get(), "element_formulation", - "B33-3D-Euler-Bernoulli"); - } - RequireStringAttribute(metadata.Get(), "step_name", kStepName); - RequireUint64Attribute(metadata.Get(), "frame_index", 0U); - - RequireCompoundDataset( - file.Get(), "/model/nodes", static_cast(domain.Nodes().size()), - {"internal_node_id", "instance_name", "source_label", "coordinates"}); - auto nodes = OpenDatasetForCheck(file.Get(), "/model/nodes"); - RequireStringAttribute(nodes.Get(), "coordinate_system", "global-cartesian"); - RequireStringAttribute(nodes.Get(), "units_label", "length"); - const std::vector nodal_dimensions = { - static_cast(domain.Nodes().size()), kDofsPerNode}; - RequireDoubleDataset(file.Get(), "/steps/Step-1/frames/0/nodal/displacement", - nodal_dimensions, "UX,UY,UZ,URX,URY,URZ", - "length,length,length,radian,radian,radian", - "global-cartesian", "nodal"); - RequireDoubleDataset( - file.Get(), "/steps/Step-1/frames/0/nodal/reaction", nodal_dimensions, - "RF1,RF2,RF3,RM1,RM2,RM3", - "force,force,force,force*length,force*length,force*length", - "global-cartesian", "nodal"); - - if (IsShellDomain(domain)) { - RequireCompoundDataset( - file.Get(), "/model/elements", - static_cast(domain.ShellElements().Size()), - {"internal_element_id", "instance_name", "source_label", - "source_element_type", "internal_formulation", "node_internal_ids", - "shell_section_internal_id", "material_internal_id"}); - auto elements = OpenDatasetForCheck(file.Get(), "/model/elements"); - RequireStringAttribute(elements.Get(), "formulation", "FESA-MITC4"); - RequireCompoundDataset( - file.Get(), "/model/shell/materials", - static_cast(domain.LinearElasticMaterials().Size()), - {"internal_material_id", "name", "E", "nu"}); - RequireCompoundDataset( - file.Get(), "/model/shell/sections", - static_cast(domain.ShellSections().Size()), - {"internal_section_id", "source_file", "source_line", "source_elset", - "material_internal_id", "thickness"}); - - std::vector directors; - std::vector 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()); - } - RequireModelDoubleDataset(file.Get(), "/model/shell/nodal_director", - {static_cast(domain.Nodes().size()), 3U}, - "D1,D2,D3", "1,1,1", "global-cartesian", "nodal", - &directors); - RequireModelDoubleDataset( - file.Get(), "/model/shell/nodal_frame", - {static_cast(domain.Nodes().size()), 3U, 3U}, "X,Y,Z", "1,1,1", - "global-cartesian", "nodal-frame", &frames); - auto nodal_frame = - OpenDatasetForCheck(file.Get(), "/model/shell/nodal_frame"); - RequireStringAttribute(nodal_frame.Get(), "axis_names", "A,B,D"); - RequireUint8Dataset(file.Get(), "/model/nodal_constraint_mask", - nodal_dimensions, model_data.constraint_mask); - RequireModelDoubleDataset(file.Get(), "/model/prescribed_displacement", - nodal_dimensions, "UX,UY,UZ,URX,URY,URZ", - "length,length,length,radian,radian,radian", - "global-cartesian", "nodal-prescribed-value", - &model_data.prescribed_displacement); - const double gauss = 1.0 / std::sqrt(3.0); - const std::vector locations{-gauss, -gauss, gauss, -gauss, - gauss, gauss, -gauss, gauss}; - RequireModelDoubleDataset(file.Get(), "/model/shell/midsurface_locations", - {4U, 2U}, "XI,ETA", "1,1", "shell-natural", - "midsurface-location", &locations); - const std::vector section_positions{-1.0, 0.0, 1.0}; - RequireModelDoubleDataset(file.Get(), "/model/shell/section_positions", - {3U}, "ZETA", "1", "shell-natural", - "section-position", §ion_positions); - auto positions = - OpenDatasetForCheck(file.Get(), "/model/shell/section_positions"); - RequireStringAttribute(positions.Get(), "position_names", - "BOTTOM,MIDDLE,TOP"); - - const hsize_t element_count = - static_cast(domain.ShellElements().Size()); - const std::string shell_root = std::string{kStepRoot} + "/element/shell"; - RequireDoubleDataset(file.Get(), (shell_root + "/local_frame").c_str(), - {element_count, 4U, 3U, 3U}, "X,Y,Z", "1,1,1", - "global-cartesian", "shell-local-frame"); - auto local_frame = - OpenDatasetForCheck(file.Get(), (shell_root + "/local_frame").c_str()); - RequireStringAttribute(local_frame.Get(), "axis_names", "E1,E2,E3"); - RequireStringAttribute(local_frame.Get(), "internal_formulation", - "FESA-MITC4"); - RequireStringAttribute(local_frame.Get(), "source_element_type_dataset", - "/model/elements.source_element_type"); - RequireStringAttribute(local_frame.Get(), "midsurface_location_dataset", - "/model/shell/midsurface_locations"); - RequireDoubleDataset( - file.Get(), (shell_root + "/generalized_strain").c_str(), - {element_count, 4U, 8U}, "E11,E22,G12,K11,K22,K12,G13,G23", - "1,1,1,1/length,1/length,1/length,1,1", "shell-local", "midsurface"); - RequireShellResultIdentity( - file.Get(), (shell_root + "/generalized_strain").c_str(), false); - RequireDoubleDataset( - file.Get(), (shell_root + "/section_resultant").c_str(), - {element_count, 4U, 8U}, "N11,N22,N12,M11,M22,M12,Q13,Q23", - "force/length,force/length,force/length,force,force,force,force/" - "length,force/length", - "shell-local", "midsurface"); - RequireShellResultIdentity( - file.Get(), (shell_root + "/section_resultant").c_str(), false); - RequireDoubleDataset(file.Get(), (shell_root + "/stress").c_str(), - {element_count, 4U, 3U, 3U}, "S11,S22,S12", - "force/length^2,force/length^2,force/length^2", - "shell-local", "section-position"); - RequireShellResultIdentity(file.Get(), (shell_root + "/stress").c_str(), - true); - RequireDoubleDataset(file.Get(), "/steps/Step-1/frames/0/global/energy", - {1U}, "PHYSICAL_STRAIN_ENERGY", "force*length", - "global", "global"); - RequireDoubleDataset( - file.Get(), "/steps/Step-1/frames/0/global/equilibrium", {6U}, - "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"); - RequireDoubleDataset( - file.Get(), "/steps/Step-1/frames/0/global/verification_metrics", {3U}, - "FREE_RESIDUAL_NORMALIZED,FORCE_BALANCE_NORMALIZED,MOMENT_BALANCE_" - "NORMALIZED", - "1,1,1", "global", "verification"); - auto metrics = OpenDatasetForCheck( - file.Get(), "/steps/Step-1/frames/0/global/verification_metrics"); - RequireStringAttribute( - metrics.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"); - RequireStringAttribute(metrics.Get(), "acceptance_thresholds", - "1e-10,1e-10,1e-10"); - for (const char* forbidden : - {"/steps/Step-1/frames/0/element/shell/drilling", - "/steps/Step-1/frames/0/element/shell/drilling_energy", - "/steps/Step-1/frames/0/element/shell/S33", - "/steps/Step-1/frames/0/element/shell/S13", - "/steps/Step-1/frames/0/element/shell/S23"}) { - if (H5Lexists(file.Get(), forbidden, H5P_DEFAULT) != 0) { - throw Hdf5Failure{"A forbidden shell result path exists."}; - } - } - } else { - RequireCompoundDataset(file.Get(), "/model/elements", - static_cast(domain.BeamElements().Size()), - {"internal_element_id", "instance_name", - "source_label", "node_internal_ids", "local_axes"}); - auto elements = OpenDatasetForCheck(file.Get(), "/model/elements"); - RequireStringAttribute(elements.Get(), "formulation", - "B33-3D-Euler-Bernoulli"); - const std::vector end_dimensions = { - static_cast(domain.BeamElements().Size()), kEndpointCount, - kEndActionComponentCount}; - const std::vector generalized_dimensions = { - static_cast(domain.BeamElements().Size()), kGaussPointCount, - kGeneralizedComponentCount}; - RequireDoubleDataset( - file.Get(), "/steps/Step-1/frames/0/element/end_force_local", - end_dimensions, "FX,FY,FZ,MX,MY,MZ", - "force,force,force,force*length,force*length,force*length", - "beam-local", "endpoint-outward-action"); - RequireDoubleDataset(file.Get(), - "/steps/Step-1/frames/0/element/section_resultant", - generalized_dimensions, "N,T,My,Mz", - "force,force*length,force*length,force*length", - "beam-local", "endpoint-positive-local-x-section-cut"); - RequireDoubleDataset( - file.Get(), "/steps/Step-1/frames/0/element/generalized_strain", - generalized_dimensions, "epsilon0,kappa_x,kappa_y,kappa_z", - "1,1/length,1/length,1/length", "beam-local", "integration-point"); - RequireDoubleDataset(file.Get(), - "/steps/Step-1/frames/0/element/generalized_resultant", - generalized_dimensions, "N,T,My,Mz", - "force,force*length,force*length,force*length", - "beam-local", "integration-point"); - RequireCompoundDataset( - file.Get(), "/steps/Step-1/frames/0/element/stress_s11", - static_cast(state.StressResults().size()), - {"internal_element_id", "gauss_point_index", "section_point_index", - "x1", "x2", "source", "S11"}); - auto stress = OpenDatasetForCheck( - file.Get(), "/steps/Step-1/frames/0/element/stress_s11"); - RequireResultAttributes(stress.Get(), "S11", "force/length^2", "beam-local", - "section-point"); - } - RequireCompoundDataset(file.Get(), "/diagnostics", - static_cast(diagnostic_count), - {"severity", "code", "file", "line", "keyword", - "entity_identity", "message"}); - RequireHdf5(file.CloseChecked(), - "Unable to close the read-only HDF5 schema self-check handle."); -} - -/// @brief Selects a unique same-directory path for atomic finalization. -std::filesystem::path MakeTemporaryPath( - const std::filesystem::path& output_path) { - static std::atomic sequence{0U}; - const std::filesystem::path parent = output_path.parent_path(); - for (std::size_t attempt = 0U; attempt < 1024U; ++attempt) { - const std::wstring filename = - L"." + output_path.filename().wstring() + L".tmp." + - std::to_wstring(GetCurrentProcessId()) + L"." + - std::to_wstring(sequence.fetch_add(1U)); - const std::filesystem::path candidate = parent / filename; - std::error_code error; - const bool exists = std::filesystem::exists(candidate, error); - if (error) { - throw Hdf5Failure{"Unable to inspect the temporary output path."}; - } - if (!exists) { - return candidate; - } - } - throw Hdf5Failure{"Unable to allocate a unique temporary output path."}; -} - -/// @brief Atomically publishes a checked candidate without weakening an -/// existing final result on failure. -bool FinalizeFile(const std::filesystem::path& temporary_path, - const std::filesystem::path& output_path) { - const DWORD attributes = GetFileAttributesW(output_path.c_str()); - if (attributes != INVALID_FILE_ATTRIBUTES) { - // ReplaceFileW is the one-operation existing-final policy: failure - // leaves the old authoritative bytes in place. - return ReplaceFileW(output_path.c_str(), temporary_path.c_str(), nullptr, - 0U, nullptr, nullptr) != FALSE; - } - const DWORD error = GetLastError(); - if (error != ERROR_FILE_NOT_FOUND && error != ERROR_PATH_NOT_FOUND) { - return false; - } - // A new final uses a same-directory write-through move after close/check. - return MoveFileExW(temporary_path.c_str(), output_path.c_str(), - MOVEFILE_WRITE_THROUGH) != FALSE; -} - -} // namespace Status Hdf5ResultsWriter::Write(const std::filesystem::path& output_path, const Domain& domain, const AnalysisState& state, const std::vector& diagnostics) { - WriterModelData model_data; - try { - const Status validation = ValidateWriterInput(output_path, domain, state, - diagnostics, model_data); - if (!validation.IsOk()) { - return validation; - } + auto plan_result = + hdf5_internal::BuildWriterPlan(output_path, domain, state, diagnostics); + if (!plan_result.HasValue()) { + return plan_result.GetStatus(); + } + hdf5_internal::WriterPlan plan = std::move(plan_result.Value()); - Hdf5ErrorSilencer silence_backend_errors; - const std::filesystem::path temporary_path = MakeTemporaryPath(output_path); - TemporaryFileGuard cleanup{temporary_path}; - WriteFile(temporary_path, domain, state, diagnostics, model_data); - // Close occurs when WriteFile returns; only the reopened read-only file - // can approve this temp artifact for authoritative replacement. - SelfCheckFile(temporary_path, domain, state, diagnostics.size(), - model_data); - if (!FinalizeFile(temporary_path, output_path)) { - return OutputFailure("hdf5-finalization-failure", - "The checked temporary HDF5 file could not " - "replace the final output."); + try { + hdf5_internal::Hdf5ErrorSilencer silence_backend_errors; + const std::filesystem::path candidate_path = + hdf5_internal::MakeCandidatePath(output_path); + hdf5_internal::CandidateFileGuard cleanup{candidate_path}; + hdf5_internal::Hdf5Handle file{ + hdf5_internal::RequireHdf5Id( + H5Fcreate(candidate_path.string().c_str(), H5F_ACC_EXCL, + H5P_DEFAULT, H5P_DEFAULT), + "Unable to create the temporary HDF5 file."), + H5Fclose}; + (void)hdf5_internal::CreateGroup(file.Get(), "/model"); + (void)hdf5_internal::CreateGroup(file.Get(), + "/steps/Step-1/frames/0/nodal"); + (void)hdf5_internal::CreateGroup(file.Get(), + "/steps/Step-1/frames/0/element"); + if (hdf5_internal::IsShellDomain(domain)) { + (void)hdf5_internal::CreateGroup(file.Get(), "/model/shell"); + (void)hdf5_internal::CreateGroup(file.Get(), + "/steps/Step-1/frames/0/element/shell"); + (void)hdf5_internal::CreateGroup(file.Get(), + "/steps/Step-1/frames/0/global"); } - cleanup.Release(); - return Status::Ok(); - } catch (const Hdf5Failure& failure) { - return OutputFailure("hdf5-write-failure", failure.what()); + hdf5_internal::WriteModel(file.Get(), domain, plan.model_data); + hdf5_internal::WriteResults(file.Get(), domain, state, diagnostics, plan); + hdf5_internal::RequireHdf5(H5Fflush(file.Get(), H5F_SCOPE_GLOBAL), + "Unable to flush the temporary HDF5 file."); + hdf5_internal::RequireHdf5( + file.CloseChecked(), + "Unable to close the temporary HDF5 file after writing."); + + const Status self_check = hdf5_internal::SelfCheck( + candidate_path, domain, state, diagnostics.size(), plan.model_data); + const Status atomic = hdf5_internal::AtomicFinalizeValidatedCandidate( + candidate_path, output_path, self_check); + if (atomic.IsOk()) { + cleanup.Release(); + } + return atomic; + } catch (const hdf5_internal::Hdf5Failure& failure) { + return hdf5_internal::OutputFailure("hdf5-write-failure", failure.what()); } catch (const std::exception& failure) { - return OutputFailure("hdf5-write-failure", failure.what()); + return hdf5_internal::OutputFailure("hdf5-write-failure", failure.what()); } } diff --git a/src/fesa/io/hdf5/hdf5_self_check.cpp b/src/fesa/io/hdf5/hdf5_self_check.cpp new file mode 100644 index 0000000..833c7aa --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_self_check.cpp @@ -0,0 +1,610 @@ +#include "io/hdf5/hdf5_self_check.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/build_info.h" +#include "io/hdf5/hdf5_raii.h" + +namespace fesa::hdf5_internal { + +Hdf5Handle OpenDatasetForCheck(const hid_t file, const char* path) { + return Hdf5Handle{RequireHdf5Id(H5Dopen2(file, path, H5P_DEFAULT), + "A mandatory HDF5 dataset is missing."), + H5Dclose}; +} + +std::vector CheckedDimensions(const hid_t dataset) { + Hdf5Handle space{RequireHdf5Id(H5Dget_space(dataset), + "Unable to inspect an HDF5 dataspace."), + H5Sclose}; + const int rank = H5Sget_simple_extent_ndims(space.Get()); + if (rank < 0) { + throw Hdf5Failure{"Unable to inspect an HDF5 dataset rank."}; + } + std::vector dimensions(static_cast(rank)); + if (rank > 0) { + RequireHdf5( + H5Sget_simple_extent_dims(space.Get(), dimensions.data(), nullptr), + "Unable to inspect HDF5 dataset dimensions."); + } + return dimensions; +} + +void RequireStringAttribute(const hid_t object, const char* name, + const std::string& expected) { + Hdf5Handle attribute{ + RequireHdf5Id(H5Aopen(object, name, H5P_DEFAULT), + "A mandatory HDF5 string attribute is missing."), + H5Aclose}; + Hdf5Handle type{RequireHdf5Id(H5Aget_type(attribute.Get()), + "Unable to inspect an HDF5 attribute type."), + H5Tclose}; + if (H5Tget_class(type.Get()) != H5T_STRING || + H5Tis_variable_str(type.Get()) <= 0 || + H5Tget_cset(type.Get()) != H5T_CSET_UTF8) { + throw Hdf5Failure{"An HDF5 string attribute is not variable-length UTF-8."}; + } + char* raw = nullptr; + RequireHdf5(H5Aread(attribute.Get(), type.Get(), &raw), + "Unable to read an HDF5 string attribute."); + const std::string actual = raw == nullptr ? std::string{} : std::string{raw}; + if (raw != nullptr) { + RequireHdf5(H5free_memory(raw), "Unable to release an HDF5 string value."); + } + if (actual != expected) { + throw Hdf5Failure{"An HDF5 string attribute has the wrong value."}; + } +} + +void RequireUint64Attribute(const hid_t object, const char* name, + const std::uint64_t expected) { + Hdf5Handle attribute{ + RequireHdf5Id(H5Aopen(object, name, H5P_DEFAULT), + "A mandatory HDF5 integer attribute is missing."), + H5Aclose}; + Hdf5Handle type{ + RequireHdf5Id(H5Aget_type(attribute.Get()), + "Unable to inspect an HDF5 integer attribute type."), + H5Tclose}; + if (H5Tget_class(type.Get()) != H5T_INTEGER || + H5Tget_size(type.Get()) != sizeof(std::uint64_t) || + H5Tget_sign(type.Get()) != H5T_SGN_NONE || + H5Tequal(type.Get(), H5T_STD_U64LE) <= 0) { + throw Hdf5Failure{"An HDF5 integer attribute is not portable uint64."}; + } + std::uint64_t actual = 0U; + RequireHdf5(H5Aread(attribute.Get(), H5T_NATIVE_UINT64, &actual), + "Unable to read an HDF5 integer attribute."); + if (actual != expected) { + throw Hdf5Failure{"An HDF5 integer attribute has the wrong value."}; + } +} + +void RequirePortableCompoundMember(const hid_t compound_type, + const unsigned index, + const std::string& name) { + Hdf5Handle member_type{ + RequireHdf5Id(H5Tget_member_type(compound_type, index), + "Unable to inspect a compound HDF5 field type."), + H5Tclose}; + const bool is_uint64 = + name == "internal_node_id" || name == "internal_element_id" || + name == "internal_material_id" || name == "internal_section_id" || + name == "shell_section_internal_id" || name == "material_internal_id" || + name == "source_line" || name == "gauss_point_index" || + name == "section_point_index" || name == "line"; + const bool is_float64 = name == "x1" || name == "x2" || name == "S11" || + name == "E" || name == "nu" || name == "thickness"; + const bool is_string = + name == "instance_name" || name == "source_label" || + name == "source_element_type" || name == "internal_formulation" || + name == "name" || name == "source_file" || name == "source_elset" || + name == "source" || name == "severity" || name == "code" || + name == "file" || name == "keyword" || name == "entity_identity" || + name == "message"; + if (is_uint64) { + if (H5Tget_class(member_type.Get()) != H5T_INTEGER || + H5Tget_size(member_type.Get()) != sizeof(std::uint64_t) || + H5Tget_sign(member_type.Get()) != H5T_SGN_NONE || + H5Tequal(member_type.Get(), H5T_STD_U64LE) <= 0) { + throw Hdf5Failure{"A compound HDF5 field is not portable uint64."}; + } + return; + } + if (is_float64) { + if (H5Tget_class(member_type.Get()) != H5T_FLOAT || + H5Tget_size(member_type.Get()) != sizeof(double) || + H5Tequal(member_type.Get(), H5T_IEEE_F64LE) <= 0) { + throw Hdf5Failure{"A compound HDF5 field is not float64."}; + } + return; + } + if (is_string) { + if (H5Tget_class(member_type.Get()) != H5T_STRING || + H5Tis_variable_str(member_type.Get()) <= 0 || + H5Tget_cset(member_type.Get()) != H5T_CSET_UTF8) { + throw Hdf5Failure{"A compound HDF5 field is not UTF-8."}; + } + return; + } + + if (H5Tget_class(member_type.Get()) != H5T_ARRAY) { + throw Hdf5Failure{"A compound HDF5 array field has the wrong type."}; + } + const int rank = H5Tget_array_ndims(member_type.Get()); + if (rank <= 0) { + throw Hdf5Failure{"Unable to inspect a compound HDF5 array rank."}; + } + std::vector dimensions(static_cast(rank)); + RequireHdf5(H5Tget_array_dims2(member_type.Get(), dimensions.data()), + "Unable to inspect compound HDF5 array dimensions."); + Hdf5Handle base_type{ + RequireHdf5Id(H5Tget_super(member_type.Get()), + "Unable to inspect a compound HDF5 array base type."), + H5Tclose}; + if (name == "node_internal_ids") { + if ((dimensions != std::vector{2U} && + dimensions != std::vector{4U}) || + H5Tequal(base_type.Get(), H5T_STD_U64LE) <= 0) { + throw Hdf5Failure{"Element connectivity is not uint64[2] or uint64[4]."}; + } + } else if (name == "coordinates") { + if (dimensions != std::vector{3U} || + H5Tequal(base_type.Get(), H5T_IEEE_F64LE) <= 0) { + throw Hdf5Failure{"Node coordinates are not float64[3]."}; + } + } else if (name == "local_axes") { + if (dimensions != std::vector{3U, 3U} || + H5Tequal(base_type.Get(), H5T_IEEE_F64LE) <= 0) { + throw Hdf5Failure{"Element local axes are not float64[3,3]."}; + } + } else { + throw Hdf5Failure{"A compound HDF5 field has an unknown type contract."}; + } +} + +void RequireResultAttributes(const hid_t dataset, const char* component_names, + const char* component_units, + const char* coordinate_system, + const char* location) { + RequireStringAttribute(dataset, "component_names", component_names); + RequireStringAttribute(dataset, "component_unit_dimensions", component_units); + RequireStringAttribute(dataset, "coordinate_system", coordinate_system); + RequireStringAttribute(dataset, "location", location); + RequireStringAttribute(dataset, "step_name", kStepName); + RequireUint64Attribute(dataset, "frame_index", 0U); +} + +void RequireShellResultIdentity(const hid_t file, const char* path, + const bool uses_section_positions) { + auto dataset = OpenDatasetForCheck(file, path); + RequireStringAttribute(dataset.Get(), "source_element_type_dataset", + "/model/elements.source_element_type"); + RequireStringAttribute(dataset.Get(), "internal_formulation", "FESA-MITC4"); + RequireStringAttribute(dataset.Get(), "midsurface_location_dataset", + "/model/shell/midsurface_locations"); + RequireStringAttribute(dataset.Get(), "local_frame_dataset", + "/steps/Step-1/frames/0/element/shell/local_frame"); + if (uses_section_positions) { + RequireStringAttribute(dataset.Get(), "section_position_dataset", + "/model/shell/section_positions"); + } +} + +void RequireDoubleDataset(const hid_t file, const char* path, + const std::vector& expected_dimensions, + const char* component_names, + const char* component_units, + const char* coordinate_system, const char* location) { + auto dataset = OpenDatasetForCheck(file, path); + if (CheckedDimensions(dataset.Get()) != expected_dimensions) { + throw Hdf5Failure{"A floating-point HDF5 dataset has the wrong shape."}; + } + Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), + "Unable to inspect an HDF5 dataset type."), + H5Tclose}; + if (H5Tget_class(type.Get()) != H5T_FLOAT || + H5Tget_size(type.Get()) != sizeof(double) || + H5Tequal(type.Get(), H5T_IEEE_F64LE) <= 0) { + throw Hdf5Failure{ + "A result dataset is not IEEE-754 float64 little-endian."}; + } + RequireResultAttributes(dataset.Get(), component_names, component_units, + coordinate_system, location); + std::size_t value_count = 1U; + for (const hsize_t dimension : expected_dimensions) { + std::size_t next = 0U; + if (!SizeProductFits(value_count, static_cast(dimension), + next)) { + throw Hdf5Failure{"A result dataset shape overflows size_t."}; + } + value_count = next; + } + std::vector values(value_count); + if (!values.empty()) { + RequireHdf5(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), + "Unable to read a result dataset during self-check."); + } + if (!std::all_of(values.begin(), values.end(), + [](const double value) { return std::isfinite(value); })) { + throw Hdf5Failure{"A result dataset contains a nonfinite value."}; + } +} + +void RequireModelDoubleDataset( + const hid_t file, const char* path, + const std::vector& expected_dimensions, + const char* component_names, const char* component_units, + const char* coordinate_system, const char* location, + const std::vector* expected_values = nullptr) { + auto dataset = OpenDatasetForCheck(file, path); + if (CheckedDimensions(dataset.Get()) != expected_dimensions) { + throw Hdf5Failure{"A model floating-point dataset has the wrong shape."}; + } + Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), + "Unable to inspect a model dataset type."), + H5Tclose}; + if (H5Tget_class(type.Get()) != H5T_FLOAT || + H5Tget_size(type.Get()) != sizeof(double) || + H5Tequal(type.Get(), H5T_IEEE_F64LE) <= 0) { + throw Hdf5Failure{"A model dataset is not IEEE-754 float64 little-endian."}; + } + RequireStringAttribute(dataset.Get(), "component_names", component_names); + RequireStringAttribute(dataset.Get(), "component_unit_dimensions", + component_units); + RequireStringAttribute(dataset.Get(), "coordinate_system", coordinate_system); + RequireStringAttribute(dataset.Get(), "location", location); + std::size_t value_count = 1U; + for (const hsize_t dimension : expected_dimensions) { + std::size_t next = 0U; + if (!SizeProductFits(value_count, static_cast(dimension), + next)) { + throw Hdf5Failure{"A model dataset shape overflows size_t."}; + } + value_count = next; + } + std::vector values(value_count); + if (!values.empty()) { + RequireHdf5(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), + "Unable to read a model dataset during self-check."); + } + if (!std::all_of(values.begin(), values.end(), + [](const double value) { return std::isfinite(value); })) { + throw Hdf5Failure{"A model dataset contains a nonfinite value."}; + } + if (expected_values != nullptr && values != *expected_values) { + throw Hdf5Failure{"A fixed model dataset has the wrong value order."}; + } +} + +void RequireUint8Dataset(const hid_t file, const char* path, + const std::vector& expected_dimensions, + const std::vector& expected_values) { + auto dataset = OpenDatasetForCheck(file, path); + if (CheckedDimensions(dataset.Get()) != expected_dimensions) { + throw Hdf5Failure{"A uint8 HDF5 dataset has the wrong shape."}; + } + Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), + "Unable to inspect a uint8 dataset type."), + H5Tclose}; + if (H5Tget_class(type.Get()) != H5T_INTEGER || + H5Tget_size(type.Get()) != sizeof(std::uint8_t) || + H5Tget_sign(type.Get()) != H5T_SGN_NONE || + H5Tequal(type.Get(), H5T_STD_U8LE) <= 0) { + throw Hdf5Failure{"A constraint mask is not portable uint8."}; + } + std::vector values(expected_values.size()); + if (!values.empty()) { + RequireHdf5(H5Dread(dataset.Get(), H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), + "Unable to read the constraint mask during self-check."); + } + if (values != expected_values) { + throw Hdf5Failure{"The constraint mask has the wrong value order."}; + } + RequireStringAttribute(dataset.Get(), "component_names", + "UX,UY,UZ,URX,URY,URZ"); + RequireStringAttribute(dataset.Get(), "value_meaning", + "0=free,1=constrained"); +} + +void RequireCompoundDataset(const hid_t file, const char* path, + const hsize_t expected_rows, + const std::vector& expected_members) { + auto dataset = OpenDatasetForCheck(file, path); + if (CheckedDimensions(dataset.Get()) != std::vector{expected_rows}) { + throw Hdf5Failure{"A compound HDF5 dataset has the wrong shape."}; + } + Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), + "Unable to inspect a compound type."), + H5Tclose}; + if (H5Tget_class(type.Get()) != H5T_COMPOUND || + H5Tget_nmembers(type.Get()) != + static_cast(expected_members.size())) { + throw Hdf5Failure{"A compound HDF5 dataset has the wrong field count."}; + } + for (std::size_t index = 0U; index < expected_members.size(); ++index) { + char* raw_name = + H5Tget_member_name(type.Get(), static_cast(index)); + if (raw_name == nullptr) { + throw Hdf5Failure{"Unable to inspect a compound HDF5 field name."}; + } + const std::string actual_name{raw_name}; + RequireHdf5(H5free_memory(raw_name), + "Unable to release a compound field name."); + if (actual_name != expected_members[index]) { + throw Hdf5Failure{"A compound HDF5 field has the wrong identity."}; + } + RequirePortableCompoundMember(type.Get(), static_cast(index), + expected_members[index]); + } +} + +/// @brief Reopens the closed candidate and verifies its required external +/// schema. +void SelfCheckFile(const std::filesystem::path& path, const Domain& domain, + const AnalysisState& state, + const std::size_t diagnostic_count, + const WriterModelData& model_data) { + if (H5Fis_hdf5(path.string().c_str()) <= 0) { + throw Hdf5Failure{"The temporary output is not an HDF5 file."}; + } + Hdf5Handle file{ + RequireHdf5Id(H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), + "Unable to reopen the temporary HDF5 file read-only."), + H5Fclose}; + Hdf5Handle metadata{ + RequireHdf5Id(H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), + "The HDF5 metadata group is missing."), + H5Gclose}; + RequireUint64Attribute(metadata.Get(), "schema_version", 0U); + RequireStringAttribute(metadata.Get(), "solver_version", + std::string{SolverVersion()}); + RequireStringAttribute(metadata.Get(), "source_input_identity", + SourceInputIdentity(domain)); + RequireStringAttribute(metadata.Get(), "unit_system_label", + "user-consistent-unspecified"); + if (IsShellDomain(domain)) { + RequireStringAttribute(metadata.Get(), "feature_id", + "linear-static-mitc4-shell"); + RequireStringAttribute( + metadata.Get(), "coordinate_convention", + "global-cartesian; shell-local=(e1,e2,e3); positive-thickness=+zeta"); + RequireStringAttribute(metadata.Get(), "internal_formulation", + "FESA-MITC4"); + RequireStringAttribute(metadata.Get(), "integration_rule", + "2x2x2-gauss; mitc4-edge-midpoint-shear"); + } else { + RequireStringAttribute(metadata.Get(), "feature_id", + "linear-static-3d-euler-beam"); + RequireStringAttribute(metadata.Get(), "coordinate_convention", + "global-cartesian; beam-local=(t,n1,t-cross-n1)"); + RequireStringAttribute(metadata.Get(), "element_formulation", + "B33-3D-Euler-Bernoulli"); + } + RequireStringAttribute(metadata.Get(), "step_name", kStepName); + RequireUint64Attribute(metadata.Get(), "frame_index", 0U); + + RequireCompoundDataset( + file.Get(), "/model/nodes", static_cast(domain.Nodes().size()), + {"internal_node_id", "instance_name", "source_label", "coordinates"}); + auto nodes = OpenDatasetForCheck(file.Get(), "/model/nodes"); + RequireStringAttribute(nodes.Get(), "coordinate_system", "global-cartesian"); + RequireStringAttribute(nodes.Get(), "units_label", "length"); + const std::vector nodal_dimensions = { + static_cast(domain.Nodes().size()), kDofsPerNode}; + RequireDoubleDataset(file.Get(), "/steps/Step-1/frames/0/nodal/displacement", + nodal_dimensions, "UX,UY,UZ,URX,URY,URZ", + "length,length,length,radian,radian,radian", + "global-cartesian", "nodal"); + RequireDoubleDataset( + file.Get(), "/steps/Step-1/frames/0/nodal/reaction", nodal_dimensions, + "RF1,RF2,RF3,RM1,RM2,RM3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "nodal"); + + if (IsShellDomain(domain)) { + RequireCompoundDataset( + file.Get(), "/model/elements", + static_cast(domain.ShellElements().Size()), + {"internal_element_id", "instance_name", "source_label", + "source_element_type", "internal_formulation", "node_internal_ids", + "shell_section_internal_id", "material_internal_id"}); + auto elements = OpenDatasetForCheck(file.Get(), "/model/elements"); + RequireStringAttribute(elements.Get(), "formulation", "FESA-MITC4"); + RequireCompoundDataset( + file.Get(), "/model/shell/materials", + static_cast(domain.LinearElasticMaterials().Size()), + {"internal_material_id", "name", "E", "nu"}); + RequireCompoundDataset( + file.Get(), "/model/shell/sections", + static_cast(domain.ShellSections().Size()), + {"internal_section_id", "source_file", "source_line", "source_elset", + "material_internal_id", "thickness"}); + + std::vector directors; + std::vector 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()); + } + RequireModelDoubleDataset(file.Get(), "/model/shell/nodal_director", + {static_cast(domain.Nodes().size()), 3U}, + "D1,D2,D3", "1,1,1", "global-cartesian", "nodal", + &directors); + RequireModelDoubleDataset( + file.Get(), "/model/shell/nodal_frame", + {static_cast(domain.Nodes().size()), 3U, 3U}, "X,Y,Z", "1,1,1", + "global-cartesian", "nodal-frame", &frames); + auto nodal_frame = + OpenDatasetForCheck(file.Get(), "/model/shell/nodal_frame"); + RequireStringAttribute(nodal_frame.Get(), "axis_names", "A,B,D"); + RequireUint8Dataset(file.Get(), "/model/nodal_constraint_mask", + nodal_dimensions, model_data.constraint_mask); + RequireModelDoubleDataset(file.Get(), "/model/prescribed_displacement", + nodal_dimensions, "UX,UY,UZ,URX,URY,URZ", + "length,length,length,radian,radian,radian", + "global-cartesian", "nodal-prescribed-value", + &model_data.prescribed_displacement); + const double gauss = 1.0 / std::sqrt(3.0); + const std::vector locations{-gauss, -gauss, gauss, -gauss, + gauss, gauss, -gauss, gauss}; + RequireModelDoubleDataset(file.Get(), "/model/shell/midsurface_locations", + {4U, 2U}, "XI,ETA", "1,1", "shell-natural", + "midsurface-location", &locations); + const std::vector section_positions{-1.0, 0.0, 1.0}; + RequireModelDoubleDataset(file.Get(), "/model/shell/section_positions", + {3U}, "ZETA", "1", "shell-natural", + "section-position", §ion_positions); + auto positions = + OpenDatasetForCheck(file.Get(), "/model/shell/section_positions"); + RequireStringAttribute(positions.Get(), "position_names", + "BOTTOM,MIDDLE,TOP"); + + const hsize_t element_count = + static_cast(domain.ShellElements().Size()); + const std::string shell_root = std::string{kStepRoot} + "/element/shell"; + RequireDoubleDataset(file.Get(), (shell_root + "/local_frame").c_str(), + {element_count, 4U, 3U, 3U}, "X,Y,Z", "1,1,1", + "global-cartesian", "shell-local-frame"); + auto local_frame = + OpenDatasetForCheck(file.Get(), (shell_root + "/local_frame").c_str()); + RequireStringAttribute(local_frame.Get(), "axis_names", "E1,E2,E3"); + RequireStringAttribute(local_frame.Get(), "internal_formulation", + "FESA-MITC4"); + RequireStringAttribute(local_frame.Get(), "source_element_type_dataset", + "/model/elements.source_element_type"); + RequireStringAttribute(local_frame.Get(), "midsurface_location_dataset", + "/model/shell/midsurface_locations"); + RequireDoubleDataset( + file.Get(), (shell_root + "/generalized_strain").c_str(), + {element_count, 4U, 8U}, "E11,E22,G12,K11,K22,K12,G13,G23", + "1,1,1,1/length,1/length,1/length,1,1", "shell-local", "midsurface"); + RequireShellResultIdentity( + file.Get(), (shell_root + "/generalized_strain").c_str(), false); + RequireDoubleDataset( + file.Get(), (shell_root + "/section_resultant").c_str(), + {element_count, 4U, 8U}, "N11,N22,N12,M11,M22,M12,Q13,Q23", + "force/length,force/length,force/length,force,force,force,force/" + "length,force/length", + "shell-local", "midsurface"); + RequireShellResultIdentity( + file.Get(), (shell_root + "/section_resultant").c_str(), false); + RequireDoubleDataset(file.Get(), (shell_root + "/stress").c_str(), + {element_count, 4U, 3U, 3U}, "S11,S22,S12", + "force/length^2,force/length^2,force/length^2", + "shell-local", "section-position"); + RequireShellResultIdentity(file.Get(), (shell_root + "/stress").c_str(), + true); + RequireDoubleDataset(file.Get(), "/steps/Step-1/frames/0/global/energy", + {1U}, "PHYSICAL_STRAIN_ENERGY", "force*length", + "global", "global"); + RequireDoubleDataset( + file.Get(), "/steps/Step-1/frames/0/global/equilibrium", {6U}, + "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"); + RequireDoubleDataset( + file.Get(), "/steps/Step-1/frames/0/global/verification_metrics", {3U}, + "FREE_RESIDUAL_NORMALIZED,FORCE_BALANCE_NORMALIZED,MOMENT_BALANCE_" + "NORMALIZED", + "1,1,1", "global", "verification"); + auto metrics = OpenDatasetForCheck( + file.Get(), "/steps/Step-1/frames/0/global/verification_metrics"); + RequireStringAttribute( + metrics.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"); + RequireStringAttribute(metrics.Get(), "acceptance_thresholds", + "1e-10,1e-10,1e-10"); + for (const char* forbidden : + {"/steps/Step-1/frames/0/element/shell/drilling", + "/steps/Step-1/frames/0/element/shell/drilling_energy", + "/steps/Step-1/frames/0/element/shell/S33", + "/steps/Step-1/frames/0/element/shell/S13", + "/steps/Step-1/frames/0/element/shell/S23"}) { + if (H5Lexists(file.Get(), forbidden, H5P_DEFAULT) != 0) { + throw Hdf5Failure{"A forbidden shell result path exists."}; + } + } + } else { + RequireCompoundDataset(file.Get(), "/model/elements", + static_cast(domain.BeamElements().Size()), + {"internal_element_id", "instance_name", + "source_label", "node_internal_ids", "local_axes"}); + auto elements = OpenDatasetForCheck(file.Get(), "/model/elements"); + RequireStringAttribute(elements.Get(), "formulation", + "B33-3D-Euler-Bernoulli"); + const std::vector end_dimensions = { + static_cast(domain.BeamElements().Size()), kEndpointCount, + kEndActionComponentCount}; + const std::vector generalized_dimensions = { + static_cast(domain.BeamElements().Size()), kGaussPointCount, + kGeneralizedComponentCount}; + RequireDoubleDataset( + file.Get(), "/steps/Step-1/frames/0/element/end_force_local", + end_dimensions, "FX,FY,FZ,MX,MY,MZ", + "force,force,force,force*length,force*length,force*length", + "beam-local", "endpoint-outward-action"); + RequireDoubleDataset(file.Get(), + "/steps/Step-1/frames/0/element/section_resultant", + generalized_dimensions, "N,T,My,Mz", + "force,force*length,force*length,force*length", + "beam-local", "endpoint-positive-local-x-section-cut"); + RequireDoubleDataset( + file.Get(), "/steps/Step-1/frames/0/element/generalized_strain", + generalized_dimensions, "epsilon0,kappa_x,kappa_y,kappa_z", + "1,1/length,1/length,1/length", "beam-local", "integration-point"); + RequireDoubleDataset(file.Get(), + "/steps/Step-1/frames/0/element/generalized_resultant", + generalized_dimensions, "N,T,My,Mz", + "force,force*length,force*length,force*length", + "beam-local", "integration-point"); + RequireCompoundDataset( + file.Get(), "/steps/Step-1/frames/0/element/stress_s11", + static_cast(state.StressResults().size()), + {"internal_element_id", "gauss_point_index", "section_point_index", + "x1", "x2", "source", "S11"}); + auto stress = OpenDatasetForCheck( + file.Get(), "/steps/Step-1/frames/0/element/stress_s11"); + RequireResultAttributes(stress.Get(), "S11", "force/length^2", "beam-local", + "section-point"); + } + RequireCompoundDataset(file.Get(), "/diagnostics", + static_cast(diagnostic_count), + {"severity", "code", "file", "line", "keyword", + "entity_identity", "message"}); + RequireHdf5(file.CloseChecked(), + "Unable to close the read-only HDF5 schema self-check handle."); +} + +Status SelfCheck(const std::filesystem::path& path, const Domain& domain, + const AnalysisState& state, const std::size_t diagnostic_count, + const WriterModelData& model_data) { + try { + SelfCheckFile(path, domain, state, diagnostic_count, model_data); + return Status::Ok(); + } catch (const Hdf5Failure& failure) { + return OutputFailure("hdf5-write-failure", failure.what()); + } catch (const std::exception& failure) { + return OutputFailure("hdf5-write-failure", failure.what()); + } +} + +} // namespace fesa::hdf5_internal diff --git a/src/fesa/io/hdf5/hdf5_self_check.h b/src/fesa/io/hdf5/hdf5_self_check.h new file mode 100644 index 0000000..b076b1e --- /dev/null +++ b/src/fesa/io/hdf5/hdf5_self_check.h @@ -0,0 +1,21 @@ +#ifndef FESA_IO_HDF5_HDF5_SELF_CHECK_H_ +#define FESA_IO_HDF5_HDF5_SELF_CHECK_H_ + +#include +#include + +#include "fesa/analysis/analysis_state.h" +#include "fesa/core/status.h" +#include "fesa/model/domain.h" +#include "io/hdf5/hdf5_primitives.h" + +namespace fesa::hdf5_internal { + +/// @brief Reopens and validates the complete closed schema-v0 candidate. +Status SelfCheck(const std::filesystem::path& path, const Domain& domain, + const AnalysisState& state, std::size_t diagnostic_count, + const WriterModelData& model_data); + +} // namespace fesa::hdf5_internal + +#endif // FESA_IO_HDF5_HDF5_SELF_CHECK_H_ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 824d5fe..8c4f4cd 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -29,6 +29,7 @@ add_executable( unit/io/abaqus/input_reader_test.cpp unit/io/abaqus/input_syntax_test.cpp unit/io/hdf5/hdf5_results_writer_test.cpp + unit/io/hdf5/hdf5_writer_components_test.cpp unit/loads/load_test.cpp unit/model/domain_test.cpp unit/model/model_types_test.cpp diff --git a/tests/unit/io/hdf5/hdf5_writer_components_test.cpp b/tests/unit/io/hdf5/hdf5_writer_components_test.cpp new file mode 100644 index 0000000..830eca3 --- /dev/null +++ b/tests/unit/io/hdf5/hdf5_writer_components_test.cpp @@ -0,0 +1,191 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/analysis/analysis_state.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/model/domain.h" +#include "io/hdf5/hdf5_atomic_file.h" +#include "io/hdf5/hdf5_result_writer.h" +#include "io/hdf5/hdf5_self_check.h" + +namespace { + +struct WriterFixture { + std::unique_ptr domain; + std::unique_ptr dofs; + std::unique_ptr state; +}; + +fesa::ModelDefinition MakeDefinition(const std::filesystem::path& source) { + fesa::ModelDefinition definition{}; + definition.source_path = source; + definition.source_content_identity = "fnv1a64:0123456789abcdef"; + definition.nodes = {{{"Beam-1", 101, "101"}, {0.0, 0.0, 0.0}, {source, 10U}}, + {{"Beam-1", 202, "202"}, {3.0, 4.0, 0.0}, {source, 11U}}}; + definition.materials = {{"Steel", 210.0e9, 0.3, {source, 20U}}}; + definition.sections = {{"General", + 0.02, + 3.0e-5, + 0.0, + 4.0e-5, + 5.0e-5, + {0.0, 0.0, 1.0}, + {{-0.1, 0.2}}, + {source, 30U}}}; + definition.elements = { + {{"Beam-1", 303, "303"}, {0U, 1U}, 0U, 0U, {source, 40U}}}; + definition.steps = {{"Step-1", {}, {}, 0.1, 1.0, 0.01, 1.0, {source, 50U}}}; + return definition; +} + +WriterFixture MakeFixture(const std::filesystem::path& source) { + auto domain_result = fesa::Domain::Create(MakeDefinition(source)); + if (!domain_result.HasValue()) { + throw std::runtime_error{"Component fixture Domain construction failed."}; + } + auto domain = + std::make_unique(std::move(domain_result.Value())); + auto model_result = fesa::AnalysisModel::Create(*domain); + if (!model_result.HasValue()) { + throw std::runtime_error{ + "Component fixture AnalysisModel construction failed."}; + } + const fesa::AnalysisModel model = std::move(model_result.Value()); + auto dofs_result = fesa::DofManager::Create(model); + if (!dofs_result.HasValue()) { + throw std::runtime_error{ + "Component fixture DofManager construction failed."}; + } + auto dofs = + std::make_unique(std::move(dofs_result.Value())); + auto state = std::make_unique( + fesa::AnalysisState::Create(*dofs, {"Step-1", 0U})); + for (std::size_t index = 0U; index < state->Displacement().Size(); ++index) { + state->Displacement()[index] = 0.25 + static_cast(index); + state->ExternalForce()[index] = 100.0 + static_cast(index); + state->InternalForce()[index] = 200.0 + static_cast(index); + state->Residual()[index] = 100.0 + static_cast(index); + state->Reaction()[index] = 100.0 + static_cast(index); + } + state->EndpointResults() = {{0U, + 0, + domain->Nodes()[0U].source_id, + {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}, + {11.0, 12.0, 13.0, 14.0}}, + {0U, + 1, + domain->Nodes()[1U].source_id, + {7.0, 8.0, 9.0, 10.0, 11.0, 12.0}, + {15.0, 16.0, 17.0, 18.0}}}; + state->GaussResults() = { + {0U, 1, {0.01, 0.02, 0.03, 0.04}, {21.0, 22.0, 23.0, 24.0}}, + {0U, 2, {0.05, 0.06, 0.07, 0.08}, {25.0, 26.0, 27.0, 28.0}}}; + state->StressResults() = {{0U, 1, 1U, -0.1, 0.2, 31.0, "input"}, + {0U, 2, 1U, -0.1, 0.2, 32.0, "input"}}; + return {std::move(domain), std::move(dofs), std::move(state)}; +} + +std::vector ReadBytes(const std::filesystem::path& path) { + std::ifstream input{path, std::ios::binary}; + return {std::istreambuf_iterator{input}, + std::istreambuf_iterator{}}; +} + +void WriteBytes(const std::filesystem::path& path, + const std::vector& bytes) { + std::ofstream output{path, std::ios::binary | std::ios::trunc}; + output.write(bytes.data(), static_cast(bytes.size())); + if (!output) { + throw std::runtime_error{"Unable to write component-test bytes."}; + } +} + +void ExpectOutputFailure(const fesa::Status& status, + const std::string& expected_code) { + ASSERT_FALSE(status.IsOk()); + EXPECT_EQ(status.Category(), fesa::FailureCategory::kOutput); + ASSERT_EQ(status.Diagnostics().size(), 1U); + EXPECT_EQ(status.Diagnostics()[0U].code, expected_code); +} + +} // namespace + +TEST(Hdf5WriterComponents, BuildsStableBackendNeutralBeamDatasetPlans) { + const auto directory = + std::filesystem::temp_directory_path() / "fesa-hdf5-component-plan"; + auto fixture = MakeFixture(directory / "beam.inp"); + + auto plan_result = fesa::hdf5_internal::BuildWriterPlan( + directory / "results.h5", *fixture.domain, *fixture.state, {}); + + ASSERT_TRUE(plan_result.HasValue()); + const auto& datasets = plan_result.Value().result_datasets; + ASSERT_EQ(datasets.size(), 6U); + EXPECT_EQ(datasets[0U].path, "/steps/Step-1/frames/0/nodal/displacement"); + EXPECT_EQ(datasets[0U].dimensions, std::vector({2U, 6U})); + EXPECT_EQ(datasets[2U].path, + "/steps/Step-1/frames/0/element/end_force_local"); + EXPECT_EQ(datasets[2U].dimensions, std::vector({1U, 2U, 6U})); + EXPECT_EQ(datasets[5U].path, + "/steps/Step-1/frames/0/element/generalized_resultant"); + EXPECT_TRUE(fesa::hdf5_internal::ValidateFiniteInventory(datasets).IsOk()); +} + +TEST(Hdf5WriterComponents, RejectsNonfiniteDatasetPlanInventory) { + fesa::hdf5_internal::DoubleDatasetPlan plan{ + "/steps/Step-1/frames/0/global/energy", + {1U}, + {std::numeric_limits::quiet_NaN()}, + "PHYSICAL_STRAIN_ENERGY", + "force*length", + "global", + "global"}; + + ExpectOutputFailure( + fesa::hdf5_internal::ValidateFiniteInventory({std::move(plan)}), + "invalid-result-dataset-plan"); +} + +TEST(Hdf5WriterComponents, + SelfCheckFailurePreservesExistingFinalThroughAtomicSeam) { + const auto directory = + std::filesystem::temp_directory_path() / "fesa-hdf5-component-atomic"; + std::error_code ignored; + std::filesystem::remove_all(directory, ignored); + ASSERT_TRUE(std::filesystem::create_directory(directory)); + const auto final = directory / "results.h5"; + const auto candidate = directory / ".results.h5.candidate"; + const std::vector sentinel{'v', 'a', 'l', 'i', 'd'}; + WriteBytes(final, sentinel); + WriteBytes(candidate, {'n', 'o', 't', '-', 'h', 'd', 'f', '5'}); + auto fixture = MakeFixture(directory / "beam.inp"); + auto plan_result = fesa::hdf5_internal::BuildWriterPlan( + final, *fixture.domain, *fixture.state, {}); + ASSERT_TRUE(plan_result.HasValue()); + + const fesa::Status self_check = + fesa::hdf5_internal::SelfCheck(candidate, *fixture.domain, *fixture.state, + 0U, plan_result.Value().model_data); + ExpectOutputFailure(self_check, "hdf5-write-failure"); + const fesa::Status finalize = + fesa::hdf5_internal::AtomicFinalizeValidatedCandidate(candidate, final, + self_check); + + ExpectOutputFailure(finalize, "hdf5-write-failure"); + EXPECT_EQ(ReadBytes(final), sentinel); + EXPECT_FALSE(std::filesystem::exists(candidate)); + std::filesystem::remove_all(directory, ignored); +}