diff --git a/include/fesa/io/hdf5/nonlinear_hdf5_results_writer.h b/include/fesa/io/hdf5/nonlinear_hdf5_results_writer.h new file mode 100644 index 0000000..691787a --- /dev/null +++ b/include/fesa/io/hdf5/nonlinear_hdf5_results_writer.h @@ -0,0 +1,22 @@ +#ifndef FESA_IO_HDF5_NONLINEAR_HDF5_RESULTS_WRITER_H_ +#define FESA_IO_HDF5_NONLINEAR_HDF5_RESULTS_WRITER_H_ + +#include "fesa/results/nonlinear_results_writer.h" + +namespace fesa { + +/// @brief Writes authoritative schema-v1 nonlinear HDF5 results atomically. +/// @note HDF5 and platform types remain private to the implementation. +class NonlinearHdf5ResultsWriter final : public NonlinearResultsWriter { + public: + /// @brief Writes and self-checks committed nonlinear frames before publish. + /// @return Success only after atomic replacement or structured output failure. + /// @note Factor-zero state remains metadata and is not an accepted frame. + Status Write(const std::filesystem::path& output_path, + const Domain& domain, const NonlinearAnalysisState& state, + const std::vector& diagnostics) override; +}; + +} // namespace fesa + +#endif // FESA_IO_HDF5_NONLINEAR_HDF5_RESULTS_WRITER_H_ diff --git a/include/fesa/results/nonlinear_results_writer.h b/include/fesa/results/nonlinear_results_writer.h new file mode 100644 index 0000000..841c16d --- /dev/null +++ b/include/fesa/results/nonlinear_results_writer.h @@ -0,0 +1,33 @@ +#ifndef FESA_RESULTS_NONLINEAR_RESULTS_WRITER_H_ +#define FESA_RESULTS_NONLINEAR_RESULTS_WRITER_H_ + +#include +#include + +#include "fesa/analysis/nonlinear_analysis_state.h" +#include "fesa/core/diagnostic.h" +#include "fesa/core/status.h" +#include "fesa/model/domain.h" + +namespace fesa { + +/// @brief Isolates nonlinear authoritative result storage from solver code. +class NonlinearResultsWriter { + public: + virtual ~NonlinearResultsWriter() = default; + + /// @brief Writes complete accepted nonlinear frames. + /// @param output_path Final authoritative path in the model input directory. + /// @param domain Immutable source and semantic model identity. + /// @param state Committed state containing every accepted increment frame. + /// @param diagnostics Deterministically ordered run diagnostics. + /// @return Success only after backend-specific finalization completes. + virtual Status Write(const std::filesystem::path& output_path, + const Domain& domain, + const NonlinearAnalysisState& state, + const std::vector& diagnostics) = 0; +}; + +} // namespace fesa + +#endif // FESA_RESULTS_NONLINEAR_RESULTS_WRITER_H_ diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index 28b617c..7a1afe6 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -33,6 +33,7 @@ add_library( io/abaqus/topology_mapper.cpp io/hdf5/hdf5_atomic_file.cpp io/hdf5/hdf5_model_writer.cpp + io/hdf5/nonlinear_hdf5_results_writer.cpp io/hdf5/hdf5_primitives.cpp io/hdf5/hdf5_result_writer.cpp io/hdf5/hdf5_results_writer.cpp diff --git a/src/fesa/io/hdf5/hdf5_primitives.cpp b/src/fesa/io/hdf5/hdf5_primitives.cpp index c58076e..1256580 100644 --- a/src/fesa/io/hdf5/hdf5_primitives.cpp +++ b/src/fesa/io/hdf5/hdf5_primitives.cpp @@ -138,6 +138,34 @@ void WriteUint64Attribute(const hid_t object, const char* name, "Unable to write an HDF5 integer attribute."); } +void WriteUint8Attribute(const hid_t object, const char* name, + const std::uint8_t value) { + Hdf5Handle space{RequireHdf5Id(H5Screate(H5S_SCALAR), + "Unable to create an attribute space."), + H5Sclose}; + Hdf5Handle attribute{ + RequireHdf5Id(H5Acreate2(object, name, H5T_STD_U8LE, space.Get(), + H5P_DEFAULT, H5P_DEFAULT), + "Unable to create an HDF5 uint8 attribute."), + H5Aclose}; + RequireHdf5(H5Awrite(attribute.Get(), H5T_NATIVE_UINT8, &value), + "Unable to write an HDF5 uint8 attribute."); +} + +void WriteDoubleAttribute(const hid_t object, const char* name, + const double value) { + Hdf5Handle space{RequireHdf5Id(H5Screate(H5S_SCALAR), + "Unable to create an attribute space."), + H5Sclose}; + Hdf5Handle attribute{ + RequireHdf5Id(H5Acreate2(object, name, H5T_IEEE_F64LE, space.Get(), + H5P_DEFAULT, H5P_DEFAULT), + "Unable to create an HDF5 float64 attribute."), + H5Aclose}; + RequireHdf5(H5Awrite(attribute.Get(), H5T_NATIVE_DOUBLE, &value), + "Unable to write an HDF5 float64 attribute."); +} + Hdf5Handle CreateGroup(const hid_t file, const char* path) { Hdf5Handle link_properties{RequireHdf5Id(H5Pcreate(H5P_LINK_CREATE), "Unable to create link properties."), diff --git a/src/fesa/io/hdf5/hdf5_primitives.h b/src/fesa/io/hdf5/hdf5_primitives.h index 429448b..da215e4 100644 --- a/src/fesa/io/hdf5/hdf5_primitives.h +++ b/src/fesa/io/hdf5/hdf5_primitives.h @@ -97,6 +97,8 @@ 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); +void WriteUint8Attribute(hid_t object, const char* name, std::uint8_t value); +void WriteDoubleAttribute(hid_t object, const char* name, double 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, diff --git a/src/fesa/io/hdf5/hdf5_self_check.cpp b/src/fesa/io/hdf5/hdf5_self_check.cpp index 833c7aa..f76ade7 100644 --- a/src/fesa/io/hdf5/hdf5_self_check.cpp +++ b/src/fesa/io/hdf5/hdf5_self_check.cpp @@ -59,7 +59,9 @@ void RequireStringAttribute(const hid_t object, const char* name, RequireHdf5(H5free_memory(raw), "Unable to release an HDF5 string value."); } if (actual != expected) { - throw Hdf5Failure{"An HDF5 string attribute has the wrong value."}; + throw Hdf5Failure{"HDF5 string attribute '" + std::string{name} + + "' mismatch: expected '" + expected + "', got '" + + actual + "'."}; } } @@ -99,9 +101,13 @@ void RequirePortableCompoundMember(const hid_t compound_type, 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"; + name == "section_point_index" || name == "line" || + name == "frame_index" || name == "increment_index" || + name == "iteration_count"; const bool is_float64 = name == "x1" || name == "x2" || name == "S11" || - name == "E" || name == "nu" || name == "thickness"; + name == "E" || name == "nu" || name == "thickness" || + name == "load_factor"; + const bool is_uint8 = name == "is_factor_one"; const bool is_string = name == "instance_name" || name == "source_label" || name == "source_element_type" || name == "internal_formulation" || @@ -126,6 +132,15 @@ void RequirePortableCompoundMember(const hid_t compound_type, } return; } + if (is_uint8) { + if (H5Tget_class(member_type.Get()) != H5T_INTEGER || + H5Tget_size(member_type.Get()) != sizeof(std::uint8_t) || + H5Tget_sign(member_type.Get()) != H5T_SGN_NONE || + H5Tequal(member_type.Get(), H5T_STD_U8LE) <= 0) { + throw Hdf5Failure{"A compound HDF5 field is not portable uint8."}; + } + return; + } if (is_string) { if (H5Tget_class(member_type.Get()) != H5T_STRING || H5Tis_variable_str(member_type.Get()) <= 0 || @@ -342,13 +357,159 @@ void RequireCompoundDataset(const hid_t file, const char* path, 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."}; + throw Hdf5Failure{"Compound HDF5 field identity mismatch: expected '" + + std::string{expected_members[index]} + "', got '" + + actual_name + "'."}; } RequirePortableCompoundMember(type.Get(), static_cast(index), expected_members[index]); } } +void RequireUint8Attribute(const hid_t object, const char* name, + const std::uint8_t expected) { + Hdf5Handle attribute{ + RequireHdf5Id(H5Aopen(object, name, H5P_DEFAULT), + "A mandatory HDF5 uint8 attribute is missing."), + H5Aclose}; + Hdf5Handle type{ + RequireHdf5Id(H5Aget_type(attribute.Get()), + "Unable to inspect an HDF5 uint8 attribute 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{"An HDF5 attribute is not portable uint8."}; + } + std::uint8_t actual = 0U; + RequireHdf5(H5Aread(attribute.Get(), H5T_NATIVE_UINT8, &actual), + "Unable to read an HDF5 uint8 attribute."); + if (actual != expected) { + throw Hdf5Failure{"An HDF5 uint8 attribute has the wrong value."}; + } +} + +void RequireDoubleAttribute(const hid_t object, const char* name, + const double expected) { + Hdf5Handle attribute{ + RequireHdf5Id(H5Aopen(object, name, H5P_DEFAULT), + "A mandatory HDF5 float64 attribute is missing."), + H5Aclose}; + Hdf5Handle type{ + RequireHdf5Id(H5Aget_type(attribute.Get()), + "Unable to inspect an HDF5 float64 attribute type."), + H5Tclose}; + if (H5Tget_class(type.Get()) != H5T_FLOAT || + H5Tget_size(type.Get()) != sizeof(double) || + H5Tequal(type.Get(), H5T_IEEE_F64LE) <= 0) { + throw Hdf5Failure{"An HDF5 attribute is not portable float64."}; + } + double actual = 0.0; + RequireHdf5(H5Aread(attribute.Get(), H5T_NATIVE_DOUBLE, &actual), + "Unable to read an HDF5 float64 attribute."); + if (actual != expected) { + throw Hdf5Failure{"An HDF5 float64 attribute has the wrong value."}; + } +} + +void RequireSchemaV1NumericDataset( + 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 AcceptedFrameState* frame) { + auto dataset = OpenDatasetForCheck(file, path); + if (CheckedDimensions(dataset.Get()) != expected_dimensions) { + throw Hdf5Failure{"A schema-v1 numeric dataset has the wrong shape."}; + } + Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), + "Unable to inspect a schema-v1 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 schema-v1 numeric dataset is not float64."}; + } + RequireUint64Attribute(dataset.Get(), "schema_version", 1U); + 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); + if (frame != nullptr) { + RequireStringAttribute(dataset.Get(), "step_name", "Step-1"); + RequireUint64Attribute(dataset.Get(), "frame_index", frame->frame_index); + RequireUint64Attribute(dataset.Get(), "increment_index", + frame->increment_index); + RequireDoubleAttribute(dataset.Get(), "load_factor", frame->load_factor); + RequireUint8Attribute(dataset.Get(), "is_factor_one", + frame->is_factor_one ? 1U : 0U); + } + 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 schema-v1 numeric 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 schema-v1 numeric dataset."); + } + if (!std::all_of(values.begin(), values.end(), + [](const double value) { return std::isfinite(value); })) { + throw Hdf5Failure{"A schema-v1 numeric dataset contains a nonfinite value."}; + } +} + +void RequireSchemaV1MaskDataset(const hid_t file, + const std::vector& dimensions) { + auto dataset = OpenDatasetForCheck(file, "/model/nodal_constraint_mask"); + if (CheckedDimensions(dataset.Get()) != dimensions) { + throw Hdf5Failure{"The schema-v1 constraint mask has the wrong shape."}; + } + Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), + "Unable to inspect the constraint mask 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{"The schema-v1 constraint mask is not uint8."}; + } + RequireUint64Attribute(dataset.Get(), "schema_version", 1U); + RequireStringAttribute(dataset.Get(), "component_names", "UX,UY,UZ,URX,URY,URZ"); + RequireStringAttribute(dataset.Get(), "component_unit_dimensions", + "1,1,1,1,1,1"); + RequireStringAttribute(dataset.Get(), "coordinate_system", "global-cartesian"); + RequireStringAttribute(dataset.Get(), "location", "nodal-constraint"); + RequireStringAttribute(dataset.Get(), "value_meaning", "0=free,1=constrained"); +} + +void RequireSchemaV1ShellIdentity(const hid_t file, const char* path, + const bool has_section_positions, + const std::string& frame_root) { + auto dataset = OpenDatasetForCheck(file, path); + RequireStringAttribute(dataset.Get(), "source_element_type_dataset", + "/model/elements.source_element_type"); + RequireStringAttribute(dataset.Get(), "internal_formulation", + "FESA-MITC4-total-lagrangian"); + RequireStringAttribute(dataset.Get(), "midsurface_location_dataset", + "/model/shell/midsurface_locations"); + RequireStringAttribute(dataset.Get(), "local_frame_dataset", + frame_root + "/element/shell/local_frame"); + if (has_section_positions) { + RequireStringAttribute(dataset.Get(), "section_position_dataset", + "/model/shell/section_positions"); + } +} + /// @brief Reopens the closed candidate and verifies its required external /// schema. void SelfCheckFile(const std::filesystem::path& path, const Domain& domain, @@ -594,6 +755,175 @@ void SelfCheckFile(const std::filesystem::path& path, const Domain& domain, "Unable to close the read-only HDF5 schema self-check handle."); } +void NonlinearSelfCheckFile(const std::filesystem::path& path, + const Domain& domain, + const NonlinearAnalysisState& state, + const std::size_t diagnostic_count) { + if (H5Fis_hdf5(path.string().c_str()) <= 0) { + throw Hdf5Failure{"The nonlinear candidate is not an HDF5 file."}; + } + Hdf5Handle file{ + RequireHdf5Id(H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), + "Unable to reopen nonlinear candidate read-only."), + H5Fclose}; + Hdf5Handle metadata{ + RequireHdf5Id(H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), + "The nonlinear metadata group is missing."), + H5Gclose}; + RequireUint64Attribute(metadata.Get(), "schema_version", 1U); + RequireStringAttribute(metadata.Get(), "feature_id", + "nonlinear-static-newton-raphson"); + 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"); + RequireStringAttribute( + metadata.Get(), "coordinate_convention", + "global-cartesian; shell-local=fixed-initial-reference; " + "positive-thickness=+zeta"); + RequireStringAttribute(metadata.Get(), "analysis_procedure", + "nonlinear-static-full-newton"); + RequireStringAttribute(metadata.Get(), "internal_formulation", + "FESA-MITC4-total-lagrangian"); + RequireStringAttribute(metadata.Get(), "integration_rule", + "2x2x2-gauss; mitc4-edge-midpoint-shear"); + RequireStringAttribute(metadata.Get(), "step_name", "Step-1"); + RequireStringAttribute(metadata.Get(), "frame_indexing", + "accepted-increment-order-zero-based"); + RequireUint8Attribute(metadata.Get(), "factor_zero_is_frame", 0U); + RequireDoubleAttribute(metadata.Get(), "factor_zero_load_factor", 0.0); + + const hsize_t node_count = static_cast(domain.Nodes().size()); + const hsize_t element_count = + static_cast(domain.ShellElements().Size()); + const std::vector nodal_dimensions{node_count, kDofsPerNode}; + RequireCompoundDataset(file.Get(), "/model/nodes", node_count, + {"internal_node_id", "instance_name", "source_label", + "coordinates"}); + RequireCompoundDataset( + file.Get(), "/model/elements", element_count, + {"internal_element_id", "instance_name", "source_label", + "source_element_type", "internal_formulation", "node_internal_ids", + "shell_section_internal_id", "material_internal_id"}); + 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"}); + RequireSchemaV1NumericDataset( + file.Get(), "/model/shell/nodal_director", {node_count, 3U}, "D1,D2,D3", + "1,1,1", "global-cartesian", "nodal-initial-director", nullptr); + RequireSchemaV1NumericDataset( + file.Get(), "/model/shell/nodal_frame", {node_count, 3U, 3U}, "X,Y,Z", + "1,1,1", "global-cartesian", "nodal-initial-frame", nullptr); + RequireSchemaV1MaskDataset(file.Get(), nodal_dimensions); + RequireSchemaV1NumericDataset( + file.Get(), "/model/prescribed_displacement_start", nodal_dimensions, + "UX,UY,UZ,URX,URY,URZ", "length,length,length,1,1,1", + "global-cartesian", "nodal-prescribed-start", nullptr); + RequireSchemaV1NumericDataset( + file.Get(), "/model/prescribed_displacement_target", nodal_dimensions, + "UX,UY,UZ,URX,URY,URZ", "length,length,length,1,1,1", + "global-cartesian", "nodal-prescribed-target", nullptr); + RequireSchemaV1NumericDataset(file.Get(), "/model/shell/midsurface_locations", + {4U, 2U}, "XI,ETA", "1,1", "shell-natural", + "midsurface-location", nullptr); + RequireSchemaV1NumericDataset(file.Get(), "/model/shell/section_positions", + {3U}, "ZETA", "1", "shell-natural", + "section-position", nullptr); + + const auto& frames = state.AcceptedFrames(); + RequireCompoundDataset(file.Get(), "/steps/Step-1/accepted_frames", + static_cast(frames.size()), + {"frame_index", "increment_index", "load_factor", + "iteration_count", "is_factor_one"}); + if (H5Lexists(file.Get(), "/steps/Step-1/trial_history", H5P_DEFAULT) != 0) { + throw Hdf5Failure{"Schema v1 must not contain nonlinear trial history."}; + } + for (const auto& frame : frames) { + const std::string root = "/steps/Step-1/frames/" + + std::to_string(frame.frame_index); + RequireSchemaV1NumericDataset( + file.Get(), (root + "/nodal/displacement").c_str(), nodal_dimensions, + "UX,UY,UZ,URX,URY,URZ", "length,length,length,1,1,1", + "global-cartesian", "nodal", &frame); + RequireSchemaV1NumericDataset( + file.Get(), (root + "/nodal/triad").c_str(), {node_count, 3U, 3U}, + "X,Y,Z", "1,1,1", "global-cartesian", "nodal-triad", &frame); + RequireSchemaV1NumericDataset( + file.Get(), (root + "/nodal/director").c_str(), {node_count, 3U}, + "D1,D2,D3", "1,1,1", "global-cartesian", "nodal-director", &frame); + RequireSchemaV1NumericDataset( + file.Get(), (root + "/nodal/generalized_internal_force").c_str(), + nodal_dimensions, "FI1,FI2,FI3,MI_PHI1,MI_PHI2,MI_PHI3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "nodal", &frame); + RequireSchemaV1NumericDataset( + file.Get(), (root + "/nodal/generalized_external_force").c_str(), + nodal_dimensions, "FE1,FE2,FE3,ME_PHI1,ME_PHI2,ME_PHI3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "nodal", &frame); + RequireSchemaV1NumericDataset( + file.Get(), (root + "/nodal/generalized_residual").c_str(), + nodal_dimensions, "RU1,RU2,RU3,RPHI1,RPHI2,RPHI3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "nodal", &frame); + RequireSchemaV1NumericDataset( + file.Get(), (root + "/nodal/reaction").c_str(), nodal_dimensions, + "RF1,RF2,RF3,RM1,RM2,RM3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "nodal", &frame); + const std::string shell_root = root + "/element/shell"; + RequireSchemaV1NumericDataset( + file.Get(), (shell_root + "/local_frame").c_str(), + {element_count, 4U, 3U, 3U}, "X,Y,Z", "1,1,1", "global-cartesian", + "shell-local-frame", &frame); + RequireSchemaV1NumericDataset( + 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", + &frame); + RequireSchemaV1ShellIdentity(file.Get(), + (shell_root + "/generalized_strain").c_str(), + false, root); + RequireSchemaV1NumericDataset( + 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", &frame); + RequireSchemaV1ShellIdentity(file.Get(), + (shell_root + "/section_resultant").c_str(), + false, root); + RequireSchemaV1NumericDataset( + 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", &frame); + RequireSchemaV1ShellIdentity(file.Get(), (shell_root + "/stress").c_str(), + true, root); + RequireSchemaV1NumericDataset(file.Get(), (root + "/global/energy").c_str(), + {1U}, "PHYSICAL_STRAIN_ENERGY", "force*length", + "global", "global", &frame); + RequireSchemaV1NumericDataset( + file.Get(), (root + "/global/equilibrium").c_str(), {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", &frame); + } + RequireCompoundDataset(file.Get(), "/diagnostics", + static_cast(diagnostic_count), + {"severity", "code", "file", "line", "keyword", + "entity_identity", "message"}); + RequireHdf5(file.CloseChecked(), + "Unable to close nonlinear read-only self-check file."); +} + Status SelfCheck(const std::filesystem::path& path, const Domain& domain, const AnalysisState& state, const std::size_t diagnostic_count, const WriterModelData& model_data) { @@ -607,4 +937,18 @@ Status SelfCheck(const std::filesystem::path& path, const Domain& domain, } } +Status NonlinearSelfCheck(const std::filesystem::path& path, + const Domain& domain, + const NonlinearAnalysisState& state, + const std::size_t diagnostic_count) { + try { + NonlinearSelfCheckFile(path, domain, state, diagnostic_count); + 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 index b076b1e..be6fb94 100644 --- a/src/fesa/io/hdf5/hdf5_self_check.h +++ b/src/fesa/io/hdf5/hdf5_self_check.h @@ -5,6 +5,7 @@ #include #include "fesa/analysis/analysis_state.h" +#include "fesa/analysis/nonlinear_analysis_state.h" #include "fesa/core/status.h" #include "fesa/model/domain.h" #include "io/hdf5/hdf5_primitives.h" @@ -16,6 +17,12 @@ Status SelfCheck(const std::filesystem::path& path, const Domain& domain, const AnalysisState& state, std::size_t diagnostic_count, const WriterModelData& model_data); +/// @brief Reopens and validates the complete closed schema-v1 candidate. +Status NonlinearSelfCheck(const std::filesystem::path& path, + const Domain& domain, + const NonlinearAnalysisState& state, + std::size_t diagnostic_count); + } // namespace fesa::hdf5_internal #endif // FESA_IO_HDF5_HDF5_SELF_CHECK_H_ diff --git a/src/fesa/io/hdf5/nonlinear_hdf5_results_writer.cpp b/src/fesa/io/hdf5/nonlinear_hdf5_results_writer.cpp new file mode 100644 index 0000000..00caafa --- /dev/null +++ b/src/fesa/io/hdf5/nonlinear_hdf5_results_writer.cpp @@ -0,0 +1,1113 @@ +#include "fesa/io/hdf5/nonlinear_hdf5_results_writer.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/build_info.h" +#include "fesa/elements/finite_rotation.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/math/vector3.h" +#include "io/hdf5/hdf5_atomic_file.h" +#include "io/hdf5/hdf5_primitives.h" +#include "io/hdf5/hdf5_raii.h" +#include "io/hdf5/hdf5_self_check.h" + +namespace fesa { +namespace { + +using hdf5_internal::Hdf5Failure; +using hdf5_internal::Hdf5Handle; + +constexpr std::size_t kDofsPerNode = 6U; +constexpr std::size_t kShellLocationsPerElement = 4U; +constexpr std::size_t kSectionPositions = 3U; + +struct NonlinearModelData { + std::vector constraint_mask; + std::vector prescribed_start; + std::vector prescribed_target; +}; + +struct NodeWriteRow { + std::uint64_t internal_node_id; + const char* instance_name; + const char* source_label; + double coordinates[3]; +}; + +struct ElementWriteRow { + std::uint64_t internal_element_id; + const char* instance_name; + const char* source_label; + const char* source_element_type; + const char* internal_formulation; + std::uint64_t node_internal_ids[4]; + std::uint64_t shell_section_internal_id; + std::uint64_t material_internal_id; +}; + +struct MaterialWriteRow { + std::uint64_t internal_material_id; + const char* name; + double youngs_modulus; + double poisson_ratio; +}; + +struct SectionWriteRow { + std::uint64_t internal_section_id; + const char* source_file; + std::uint64_t source_line; + const char* source_elset; + std::uint64_t material_internal_id; + double thickness; +}; + +struct AcceptedFrameWriteRow { + std::uint64_t frame_index; + std::uint64_t increment_index; + double load_factor; + std::uint64_t iteration_count; + std::uint8_t is_factor_one; +}; + +struct DiagnosticWriteRow { + const char* severity; + const char* code; + const char* file; + std::uint64_t line; + const char* keyword; + const char* entity_identity; + const char* message; +}; + +template +bool IsFinite(const std::array& values) { + return std::all_of(values.begin(), values.end(), + [](const double value) { return std::isfinite(value); }); +} + +bool IsFinite(const Vector& values) { + for (std::size_t index = 0U; index < values.Size(); ++index) { + if (!std::isfinite(values[index])) { + return false; + } + } + return true; +} + +bool SameIdentity(const SourceEntityId& left, const SourceEntityId& right) { + return left.instance_name == right.instance_name && + left.source_label == right.source_label && + left.source_label_text == right.source_label_text; +} + +bool IsFinite(const ShellResultRow& row) { + if (!IsFinite(row.natural_coordinates) || !IsFinite(row.generalized_strain) || + !IsFinite(row.section_resultant)) { + return false; + } + for (const auto& axis : row.local_frame) { + if (!IsFinite(axis)) { + return false; + } + } + for (const auto& stress : row.stress) { + if (!std::isfinite(stress.zeta) || !IsFinite(stress.components)) { + return false; + } + } + return true; +} + +Status Failure(const char* code, const std::string& message) { + return hdf5_internal::OutputFailure(code, message); +} + +Status ValidateDomain(const Domain& domain) { + if (!domain.BeamElements().Empty() || domain.ShellElements().Empty() || + domain.Steps().Size() != 1U || !domain.Steps()[0U].IsNonlinearStatic() || + domain.Nodes().empty() || domain.SourcePath().empty() || + domain.SourceContentIdentity().empty() || + !hdf5_internal::IsValidUtf8(domain.SourceContentIdentity())) { + return Failure("invalid-result-identity", + "Schema v1 requires one nonlinear shell Domain identity."); + } + if (domain.ShellNodeInitialFrames().size() != domain.Nodes().size()) { + return Failure("invalid-result-identity", + "Every nonlinear output node requires an initial frame."); + } + for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { + const Node& source = domain.Nodes()[node]; + const auto& frame = domain.ShellNodeInitialFrames()[node]; + if (source.source_id.source_label <= 0 || + source.source_id.source_label_text.empty() || + !hdf5_internal::IsValidUtf8(source.source_id.instance_name) || + !hdf5_internal::IsValidUtf8(source.source_id.source_label_text) || + !Vector3{source.coordinates}.IsFinite() || + frame.node_index != static_cast(node) || + !IsFinite(frame.tangent_a) || !IsFinite(frame.tangent_b) || + !IsFinite(frame.director)) { + return Failure("invalid-result-identity", + "Nonlinear node identity or initial frame is invalid."); + } + } + for (std::size_t element = 0U; element < domain.ShellElements().Size(); + ++element) { + const auto& source = domain.ShellElements()[element]; + if (source.source_id.source_label <= 0 || + source.source_id.source_label_text.empty() || + !hdf5_internal::IsValidUtf8(source.source_id.instance_name) || + !hdf5_internal::IsValidUtf8(source.source_id.source_label_text) || + source.node_indices.size() != 4U || + source.section_index >= domain.ShellSections().Size() || + source.material_index >= domain.LinearElasticMaterials().Size() || + domain.ShellSections()[source.section_index].material_index != + source.material_index) { + return Failure("invalid-result-identity", + "Nonlinear shell element identity is invalid."); + } + for (const EntityIndex node : source.node_indices) { + if (node >= domain.Nodes().size()) { + return Failure("invalid-result-identity", + "Nonlinear shell connectivity is out of range."); + } + } + } + return Status::Ok(); +} + +Status BuildModelData(const Domain& domain, NonlinearModelData& model_data) { + auto model_result = AnalysisModel::Create(domain); + if (!model_result.HasValue()) { + return Failure("invalid-result-state", + "The nonlinear writer could not create an AnalysisModel."); + } + auto dofs_result = DofManager::Create(model_result.Value()); + if (!dofs_result.HasValue()) { + return Failure("invalid-result-state", + "The nonlinear writer could not create a DofManager."); + } + const DofManager dofs = std::move(dofs_result.Value()); + std::size_t full_dof_count = 0U; + if (!hdf5_internal::SizeProductFits(domain.Nodes().size(), kDofsPerNode, + full_dof_count) || + dofs.FullDofCount() != full_dof_count || + dofs.ConstrainedDofs().size() != dofs.PrescribedValues().Size()) { + return Failure("invalid-result-state", + "Nonlinear constraint identities have an invalid size."); + } + model_data.constraint_mask.assign(full_dof_count, 0U); + model_data.prescribed_start.assign(full_dof_count, 0.0); + model_data.prescribed_target.assign(full_dof_count, 0.0); + for (std::size_t index = 0U; index < dofs.ConstrainedDofs().size(); ++index) { + const std::size_t full_dof = dofs.ConstrainedDofs()[index]; + const double target = dofs.PrescribedValues()[index]; + if (full_dof >= full_dof_count || !std::isfinite(target)) { + return Failure("invalid-result-state", + "Nonlinear prescribed values must be finite and in range."); + } + model_data.constraint_mask[full_dof] = 1U; + model_data.prescribed_target[full_dof] = target; + } + return Status::Ok(); +} + +Status ValidateState(const Domain& domain, const NonlinearAnalysisState& state) { + const auto& frames = state.AcceptedFrames(); + if (frames.empty()) { + return Failure("invalid-result-state", + "Schema v1 requires at least one accepted frame."); + } + std::size_t full_dof_count = 0U; + std::size_t shell_row_count = 0U; + if (!hdf5_internal::SizeProductFits(domain.Nodes().size(), kDofsPerNode, + full_dof_count) || + !hdf5_internal::SizeProductFits(domain.ShellElements().Size(), + kShellLocationsPerElement, + shell_row_count)) { + return Failure("invalid-result-state", + "Nonlinear result dimensions overflow size_t."); + } + const double gauss = 1.0 / std::sqrt(3.0); + const std::array + locations{ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2, + ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4}; + const std::array, kShellLocationsPerElement> + coordinates{{{-gauss, -gauss}, {gauss, -gauss}, + {gauss, gauss}, {-gauss, gauss}}}; + for (std::size_t frame_index = 0U; frame_index < frames.size(); + ++frame_index) { + const AcceptedFrameState& frame = frames[frame_index]; + const bool final_frame = frame_index + 1U == frames.size(); + if (frame.frame_index != frame_index || frame.increment_index == 0U || + !std::isfinite(frame.load_factor) || frame.load_factor <= 0.0 || + frame.load_factor > 1.0 || frame.iteration_count == 0U || + frame.is_factor_one != final_frame || + (final_frame && frame.load_factor != 1.0) || + (!final_frame && frame.load_factor >= 1.0) || + (frame_index > 0U && + (frame.increment_index <= frames[frame_index - 1U].increment_index || + frame.load_factor <= frames[frame_index - 1U].load_factor))) { + return Failure("invalid-result-state", + "Accepted-frame ordering or factor-one identity is invalid."); + } + if (frame.nodal_internal_order.size() != domain.Nodes().size() || + frame.nodal_source_ids.size() != domain.Nodes().size() || + frame.nodal_triad.size() != domain.Nodes().size() || + frame.nodal_director.size() != domain.Nodes().size() || + frame.shell_results.size() != shell_row_count || + !std::isfinite(frame.physical_strain_energy) || + !IsFinite(frame.global_equilibrium)) { + return Failure("invalid-result-state", + "Accepted-frame identity or physical inventory is invalid."); + } + const std::array vectors{ + &frame.generalized_coordinates, &frame.displacement, + &frame.rotation_vector, &frame.internal_force, + &frame.external_force, &frame.residual, + &frame.reaction}; + for (const Vector* values : vectors) { + if (values->Size() != full_dof_count || !IsFinite(*values)) { + return Failure("invalid-result-state", + "Accepted-frame vectors must be finite node-by-six data."); + } + } + for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { + if (frame.nodal_internal_order[node] != node || + !SameIdentity(frame.nodal_source_ids[node], + domain.Nodes()[node].source_id) || + !IsFinite(frame.nodal_triad[node][0U]) || + !IsFinite(frame.nodal_triad[node][1U]) || + !IsFinite(frame.nodal_triad[node][2U]) || + !IsFinite(frame.nodal_director[node])) { + return Failure("invalid-result-state", + "Accepted-frame nodal identity is invalid."); + } + } + for (std::size_t row_index = 0U; row_index < shell_row_count; + ++row_index) { + const ShellResultRow& row = frame.shell_results[row_index]; + const std::size_t element = row_index / kShellLocationsPerElement; + const std::size_t location = row_index % kShellLocationsPerElement; + if (row.element != element || row.location != locations[location] || + row.natural_coordinates != coordinates[location] || !IsFinite(row)) { + return Failure("invalid-result-state", + "Accepted shell rows must preserve element and GP order."); + } + } + } + return Status::Ok(); +} + +void WriteCommonNumericAttributes(const hid_t dataset, + const std::string& component_names, + const std::string& component_units, + const std::string& coordinate_system, + const std::string& location) { + hdf5_internal::WriteUint64Attribute(dataset, "schema_version", 1U); + hdf5_internal::WriteStringAttribute(dataset, "component_names", + component_names); + hdf5_internal::WriteStringAttribute(dataset, "component_unit_dimensions", + component_units); + hdf5_internal::WriteStringAttribute(dataset, "coordinate_system", + coordinate_system); + hdf5_internal::WriteStringAttribute(dataset, "location", location); +} + +void WriteFrameAttributes(const hid_t dataset, const AcceptedFrameState& frame) { + hdf5_internal::WriteStringAttribute(dataset, "step_name", "Step-1"); + hdf5_internal::WriteUint64Attribute(dataset, "frame_index", frame.frame_index); + hdf5_internal::WriteUint64Attribute(dataset, "increment_index", + frame.increment_index); + hdf5_internal::WriteDoubleAttribute(dataset, "load_factor", + frame.load_factor); + hdf5_internal::WriteUint8Attribute(dataset, "is_factor_one", + frame.is_factor_one ? 1U : 0U); +} + +Hdf5Handle WriteModelDoubleDataset( + const hid_t file, const std::string& path, + const std::vector& dimensions, const double* values, + const std::size_t value_count, const std::string& component_names, + const std::string& component_units, const std::string& coordinate_system, + const std::string& location) { + auto dataset = hdf5_internal::WriteDoubleValues(file, path, dimensions, values, + value_count); + WriteCommonNumericAttributes(dataset.Get(), component_names, component_units, + coordinate_system, location); + return dataset; +} + +Hdf5Handle WriteFrameDoubleDataset( + const hid_t file, const std::string& path, + const std::vector& dimensions, const double* values, + const std::size_t value_count, const std::string& component_names, + const std::string& component_units, const std::string& coordinate_system, + const std::string& location, const AcceptedFrameState& frame) { + auto dataset = WriteModelDoubleDataset(file, path, dimensions, values, + value_count, component_names, + component_units, coordinate_system, + location); + WriteFrameAttributes(dataset.Get(), frame); + return dataset; +} + +void WriteMetadata(const hid_t file, const Domain& domain) { + auto metadata = hdf5_internal::CreateGroup(file, "/metadata"); + hdf5_internal::WriteUint64Attribute(metadata.Get(), "schema_version", 1U); + hdf5_internal::WriteStringAttribute(metadata.Get(), "feature_id", + "nonlinear-static-newton-raphson"); + hdf5_internal::WriteStringAttribute(metadata.Get(), "solver_version", + std::string{SolverVersion()}); + hdf5_internal::WriteStringAttribute(metadata.Get(), "source_input_identity", + hdf5_internal::SourceInputIdentity(domain)); + hdf5_internal::WriteStringAttribute(metadata.Get(), "unit_system_label", + "user-consistent-unspecified"); + hdf5_internal::WriteStringAttribute( + metadata.Get(), "coordinate_convention", + "global-cartesian; shell-local=fixed-initial-reference; " + "positive-thickness=+zeta"); + hdf5_internal::WriteStringAttribute(metadata.Get(), "analysis_procedure", + "nonlinear-static-full-newton"); + hdf5_internal::WriteStringAttribute(metadata.Get(), "internal_formulation", + "FESA-MITC4-total-lagrangian"); + hdf5_internal::WriteStringAttribute(metadata.Get(), "integration_rule", + "2x2x2-gauss; mitc4-edge-midpoint-shear"); + hdf5_internal::WriteStringAttribute(metadata.Get(), "step_name", "Step-1"); + hdf5_internal::WriteStringAttribute(metadata.Get(), "frame_indexing", + "accepted-increment-order-zero-based"); + hdf5_internal::WriteUint8Attribute(metadata.Get(), "factor_zero_is_frame", + 0U); + hdf5_internal::WriteDoubleAttribute(metadata.Get(), "factor_zero_load_factor", + 0.0); +} + +void WriteNodes(const hid_t file, const Domain& domain) { + std::vector 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 = hdf5_internal::MakeUtf8StringType(); + const hsize_t coordinate_dimensions[] = {3U}; + Hdf5Handle file_coordinates{hdf5_internal::RequireHdf5Id( + H5Tarray_create2(H5T_IEEE_F64LE, 1, + coordinate_dimensions), + "Unable to create node coordinate type."), + H5Tclose}; + Hdf5Handle memory_coordinates{hdf5_internal::RequireHdf5Id( + H5Tarray_create2(H5T_NATIVE_DOUBLE, 1, + coordinate_dimensions), + "Unable to create node memory type."), + H5Tclose}; + Hdf5Handle file_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)), + "Unable to create nonlinear node file type."), + H5Tclose}; + Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)), + "Unable to create nonlinear node memory type."), + H5Tclose}; + const auto insert = [&](const hid_t type, const hid_t integer_type, + const hid_t coordinate_type) { + hdf5_internal::RequireHdf5( + H5Tinsert(type, "internal_node_id", HOFFSET(NodeWriteRow, internal_node_id), + integer_type), + "Unable to define nonlinear node ID."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "instance_name", HOFFSET(NodeWriteRow, instance_name), + string_type.Get()), + "Unable to define nonlinear node instance."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "source_label", HOFFSET(NodeWriteRow, source_label), + string_type.Get()), + "Unable to define nonlinear node source label."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "coordinates", HOFFSET(NodeWriteRow, coordinates), + coordinate_type), + "Unable to define nonlinear node coordinates."); + }; + insert(file_type.Get(), H5T_STD_U64LE, file_coordinates.Get()); + insert(memory_type.Get(), H5T_NATIVE_UINT64, memory_coordinates.Get()); + auto dataset = hdf5_internal::WriteCompoundDataset( + file, "/model/nodes", rows.size(), file_type.Get(), memory_type.Get(), + rows.data()); + hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U); + hdf5_internal::WriteStringAttribute(dataset.Get(), "coordinate_system", + "global-cartesian"); + hdf5_internal::WriteStringAttribute(dataset.Get(), "component_names", "X,Y,Z"); + hdf5_internal::WriteStringAttribute(dataset.Get(), "component_unit_dimensions", + "length,length,length"); + hdf5_internal::WriteStringAttribute(dataset.Get(), "location", + "nodal-reference"); +} + +void WriteElements(const hid_t file, const Domain& domain) { + std::vector 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(), + hdf5_internal::ShellSourceTypeName(element.source_type), + "FESA-MITC4-total-lagrangian", + {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 = hdf5_internal::MakeUtf8StringType(); + const hsize_t nodes[] = {4U}; + Hdf5Handle file_nodes{hdf5_internal::RequireHdf5Id( + H5Tarray_create2(H5T_STD_U64LE, 1, nodes), + "Unable to create shell connectivity type."), + H5Tclose}; + Hdf5Handle memory_nodes{hdf5_internal::RequireHdf5Id( + H5Tarray_create2(H5T_NATIVE_UINT64, 1, nodes), + "Unable to create shell connectivity memory type."), + H5Tclose}; + Hdf5Handle file_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)), + "Unable to create shell element file type."), + H5Tclose}; + Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)), + "Unable to create shell element memory type."), + H5Tclose}; + const auto insert = [&](const hid_t type, const hid_t integer_type, + const hid_t node_type) { + hdf5_internal::RequireHdf5( + H5Tinsert(type, "internal_element_id", + HOFFSET(ElementWriteRow, internal_element_id), integer_type), + "Unable to define shell element ID."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "instance_name", HOFFSET(ElementWriteRow, instance_name), + string_type.Get()), + "Unable to define shell instance identity."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "source_label", HOFFSET(ElementWriteRow, source_label), + string_type.Get()), + "Unable to define shell source label identity."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "source_element_type", + HOFFSET(ElementWriteRow, source_element_type), + string_type.Get()), + "Unable to define shell source type identity."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "internal_formulation", + HOFFSET(ElementWriteRow, internal_formulation), + string_type.Get()), + "Unable to define shell formulation identity."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "node_internal_ids", HOFFSET(ElementWriteRow, node_internal_ids), + node_type), + "Unable to define shell connectivity."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "shell_section_internal_id", + HOFFSET(ElementWriteRow, shell_section_internal_id), integer_type), + "Unable to define shell section identity."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "material_internal_id", + HOFFSET(ElementWriteRow, material_internal_id), integer_type), + "Unable to define shell material identity."); + }; + insert(file_type.Get(), H5T_STD_U64LE, file_nodes.Get()); + insert(memory_type.Get(), H5T_NATIVE_UINT64, memory_nodes.Get()); + auto dataset = hdf5_internal::WriteCompoundDataset( + file, "/model/elements", rows.size(), file_type.Get(), memory_type.Get(), + rows.data()); + hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U); + hdf5_internal::WriteStringAttribute(dataset.Get(), "formulation", + "FESA-MITC4-total-lagrangian"); +} + +void WriteMaterials(const hid_t file, const Domain& domain) { + std::vector 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 = hdf5_internal::MakeUtf8StringType(); + Hdf5Handle file_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(MaterialWriteRow)), + "Unable to create shell material file type."), + H5Tclose}; + Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(MaterialWriteRow)), + "Unable to create shell material memory type."), + H5Tclose}; + const auto insert = [&](const hid_t type, const hid_t integer_type, + const hid_t double_type) { + hdf5_internal::RequireHdf5( + H5Tinsert(type, "internal_material_id", + HOFFSET(MaterialWriteRow, internal_material_id), integer_type), + "Unable to define shell material ID."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "name", HOFFSET(MaterialWriteRow, name), string_type.Get()), + "Unable to define shell material name."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "E", HOFFSET(MaterialWriteRow, youngs_modulus), + double_type), + "Unable to define shell material E."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "nu", HOFFSET(MaterialWriteRow, poisson_ratio), + double_type), + "Unable to define shell material nu."); + }; + insert(file_type.Get(), H5T_STD_U64LE, H5T_IEEE_F64LE); + insert(memory_type.Get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE); + auto dataset = hdf5_internal::WriteCompoundDataset( + file, "/model/shell/materials", rows.size(), file_type.Get(), + memory_type.Get(), rows.data()); + hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U); + hdf5_internal::WriteStringAttribute(dataset.Get(), "component_unit_dimensions", + "force/length^2,1"); + hdf5_internal::WriteStringAttribute(dataset.Get(), "location", + "model-material"); +} + +void WriteSections(const hid_t file, const Domain& domain) { + std::vector source_files; + source_files.reserve(domain.ShellSections().Size()); + for (std::size_t index = 0U; index < domain.ShellSections().Size(); ++index) { + source_files.push_back(hdf5_internal::NormalizedPathString( + domain.ShellSections()[index].location.file)); + } + std::vector 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 = hdf5_internal::MakeUtf8StringType(); + Hdf5Handle file_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(SectionWriteRow)), + "Unable to create shell section file type."), + H5Tclose}; + Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(SectionWriteRow)), + "Unable to create shell section memory type."), + H5Tclose}; + const auto insert = [&](const hid_t type, const hid_t integer_type, + const hid_t double_type) { + hdf5_internal::RequireHdf5( + H5Tinsert(type, "internal_section_id", + HOFFSET(SectionWriteRow, internal_section_id), integer_type), + "Unable to define shell section ID."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "source_file", HOFFSET(SectionWriteRow, source_file), + string_type.Get()), + "Unable to define shell section source file."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "source_line", HOFFSET(SectionWriteRow, source_line), + integer_type), + "Unable to define shell section source line."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "source_elset", HOFFSET(SectionWriteRow, source_elset), + string_type.Get()), + "Unable to define shell section source elset."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "material_internal_id", + HOFFSET(SectionWriteRow, material_internal_id), integer_type), + "Unable to define shell section material identity."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "thickness", HOFFSET(SectionWriteRow, thickness), + double_type), + "Unable to define shell section thickness."); + }; + insert(file_type.Get(), H5T_STD_U64LE, H5T_IEEE_F64LE); + insert(memory_type.Get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE); + auto dataset = hdf5_internal::WriteCompoundDataset( + file, "/model/shell/sections", rows.size(), file_type.Get(), + memory_type.Get(), rows.data()); + hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U); + hdf5_internal::WriteStringAttribute(dataset.Get(), "component_unit_dimensions", + "length"); + hdf5_internal::WriteStringAttribute(dataset.Get(), "location", + "model-shell-section"); + hdf5_internal::WriteStringAttribute(dataset.Get(), "layering", + "centered-single-layer"); +} + +void WriteModelData(const hid_t file, const Domain& domain, + const NonlinearModelData& model_data) { + std::vector 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()); + } + const hsize_t node_count = static_cast(domain.Nodes().size()); + (void)WriteModelDoubleDataset(file, "/model/shell/nodal_director", + {node_count, 3U}, directors.data(), + directors.size(), "D1,D2,D3", "1,1,1", + "global-cartesian", "nodal-initial-director"); + auto initial_frame = WriteModelDoubleDataset( + file, "/model/shell/nodal_frame", {node_count, 3U, 3U}, frames.data(), + frames.size(), "X,Y,Z", "1,1,1", "global-cartesian", + "nodal-initial-frame"); + hdf5_internal::WriteStringAttribute(initial_frame.Get(), "axis_names", "A,B,D"); + + WriteMaterials(file, domain); + WriteSections(file, domain); + const std::vector nodal_dimensions{node_count, kDofsPerNode}; + hdf5_internal::WriteUint8Dataset(file, "/model/nodal_constraint_mask", + nodal_dimensions, model_data.constraint_mask.data(), + model_data.constraint_mask.size()); + Hdf5Handle mask{hdf5_internal::RequireHdf5Id( + H5Dopen2(file, "/model/nodal_constraint_mask", H5P_DEFAULT), + "Unable to reopen nonlinear constraint mask."), + H5Dclose}; + hdf5_internal::WriteUint64Attribute(mask.Get(), "schema_version", 1U); + hdf5_internal::WriteStringAttribute(mask.Get(), "component_unit_dimensions", + "1,1,1,1,1,1"); + hdf5_internal::WriteStringAttribute(mask.Get(), "coordinate_system", + "global-cartesian"); + hdf5_internal::WriteStringAttribute(mask.Get(), "location", + "nodal-constraint"); + auto start = WriteModelDoubleDataset( + file, "/model/prescribed_displacement_start", nodal_dimensions, + model_data.prescribed_start.data(), model_data.prescribed_start.size(), + "UX,UY,UZ,URX,URY,URZ", "length,length,length,1,1,1", + "global-cartesian", "nodal-prescribed-start"); + hdf5_internal::WriteDoubleAttribute(start.Get(), "load_factor", 0.0); + auto target = WriteModelDoubleDataset( + file, "/model/prescribed_displacement_target", nodal_dimensions, + model_data.prescribed_target.data(), model_data.prescribed_target.size(), + "UX,UY,UZ,URX,URY,URZ", "length,length,length,1,1,1", + "global-cartesian", "nodal-prescribed-target"); + hdf5_internal::WriteDoubleAttribute(target.Get(), "load_factor", 1.0); + + const double gauss = 1.0 / std::sqrt(3.0); + const std::array locations{-gauss, -gauss, gauss, -gauss, + gauss, gauss, -gauss, gauss}; + (void)WriteModelDoubleDataset(file, "/model/shell/midsurface_locations", + {4U, 2U}, locations.data(), locations.size(), + "XI,ETA", "1,1", "shell-natural", + "midsurface-location"); + const std::array positions{-1.0, 0.0, 1.0}; + auto section_positions = WriteModelDoubleDataset( + file, "/model/shell/section_positions", {3U}, positions.data(), + positions.size(), "ZETA", "1", "shell-natural", "section-position"); + hdf5_internal::WriteStringAttribute(section_positions.Get(), "position_names", + "BOTTOM,MIDDLE,TOP"); +} + +void WriteAcceptedFrameTable(const hid_t file, + const std::vector& frames) { + std::vector rows; + rows.reserve(frames.size()); + for (const auto& frame : frames) { + rows.push_back({frame.frame_index, frame.increment_index, frame.load_factor, + frame.iteration_count, + frame.is_factor_one ? static_cast(1U) + : static_cast(0U)}); + } + Hdf5Handle file_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(AcceptedFrameWriteRow)), + "Unable to create accepted-frame file type."), + H5Tclose}; + Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(AcceptedFrameWriteRow)), + "Unable to create accepted-frame memory type."), + H5Tclose}; + const auto insert = [&](const hid_t type, const hid_t integer_type, + const hid_t double_type, const hid_t flag_type) { + hdf5_internal::RequireHdf5( + H5Tinsert(type, "frame_index", HOFFSET(AcceptedFrameWriteRow, frame_index), + integer_type), + "Unable to define accepted-frame index."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "increment_index", + HOFFSET(AcceptedFrameWriteRow, increment_index), integer_type), + "Unable to define accepted increment index."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "load_factor", HOFFSET(AcceptedFrameWriteRow, load_factor), + double_type), + "Unable to define accepted-frame load factor."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "iteration_count", + HOFFSET(AcceptedFrameWriteRow, iteration_count), integer_type), + "Unable to define accepted-frame iteration count."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "is_factor_one", + HOFFSET(AcceptedFrameWriteRow, is_factor_one), flag_type), + "Unable to define accepted-frame factor-one flag."); + }; + insert(file_type.Get(), H5T_STD_U64LE, H5T_IEEE_F64LE, H5T_STD_U8LE); + insert(memory_type.Get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE, + H5T_NATIVE_UINT8); + auto dataset = hdf5_internal::WriteCompoundDataset( + file, "/steps/Step-1/accepted_frames", rows.size(), file_type.Get(), + memory_type.Get(), rows.data()); + hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U); +} + +std::vector 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 NodalDisplacement(const AcceptedFrameState& frame) { + std::vector values(frame.displacement.Size()); + for (std::size_t node = 0U; node < frame.nodal_internal_order.size(); ++node) { + const std::size_t offset = kDofsPerNode * node; + for (std::size_t component = 0U; component < 3U; ++component) { + values[offset + component] = frame.displacement[offset + component]; + values[offset + 3U + component] = + frame.rotation_vector[offset + 3U + component]; + } + } + return values; +} + +std::vector FlattenTriads(const AcceptedFrameState& frame) { + std::vector values; + values.reserve(frame.nodal_triad.size() * 9U); + for (const auto& triad : frame.nodal_triad) { + for (const auto& axis : triad) { + values.insert(values.end(), axis.begin(), axis.end()); + } + } + return values; +} + +std::vector FlattenDirectors(const AcceptedFrameState& frame) { + std::vector values; + values.reserve(frame.nodal_director.size() * 3U); + for (const auto& director : frame.nodal_director) { + values.insert(values.end(), director.begin(), director.end()); + } + return values; +} + +std::vector SpatialReaction(const AcceptedFrameState& frame) { + std::vector values = CopyVector(frame.reaction); + for (std::size_t node = 0U; node < frame.nodal_internal_order.size(); ++node) { + const std::size_t offset = kDofsPerNode * node; + const Vector3 phi{frame.rotation_vector[offset + 3U], + frame.rotation_vector[offset + 4U], + frame.rotation_vector[offset + 5U]}; + const Matrix transform = InverseTransposeLeftJacobianSO3(phi); + const std::array generalized{ + frame.reaction[offset + 3U], frame.reaction[offset + 4U], + frame.reaction[offset + 5U]}; + for (std::size_t row = 0U; row < 3U; ++row) { + values[offset + 3U + row] = 0.0; + for (std::size_t column = 0U; column < 3U; ++column) { + values[offset + 3U + row] += transform(row, column) * generalized[column]; + } + } + } + return values; +} + +void WriteShellIdentity(const hid_t dataset, const bool has_section_positions, + const std::string& local_frame_path) { + hdf5_internal::WriteStringAttribute(dataset, "source_element_type_dataset", + "/model/elements.source_element_type"); + hdf5_internal::WriteStringAttribute(dataset, "internal_formulation", + "FESA-MITC4-total-lagrangian"); + hdf5_internal::WriteStringAttribute(dataset, "midsurface_location_dataset", + "/model/shell/midsurface_locations"); + hdf5_internal::WriteStringAttribute(dataset, "local_frame_dataset", + local_frame_path); + if (has_section_positions) { + hdf5_internal::WriteStringAttribute(dataset, "section_position_dataset", + "/model/shell/section_positions"); + } +} + +void WriteFrame(const hid_t file, const Domain& domain, + const AcceptedFrameState& frame) { + const std::string root = "/steps/Step-1/frames/" + + std::to_string(frame.frame_index); + (void)hdf5_internal::CreateGroup(file, (root + "/nodal").c_str()); + (void)hdf5_internal::CreateGroup(file, (root + "/element/shell").c_str()); + (void)hdf5_internal::CreateGroup(file, (root + "/global").c_str()); + const hsize_t node_count = static_cast(domain.Nodes().size()); + const std::vector nodal_dimensions{node_count, kDofsPerNode}; + const auto displacement = NodalDisplacement(frame); + (void)WriteFrameDoubleDataset( + file, root + "/nodal/displacement", nodal_dimensions, displacement.data(), + displacement.size(), "UX,UY,UZ,URX,URY,URZ", "length,length,length,1,1,1", + "global-cartesian", "nodal", frame); + const auto triads = FlattenTriads(frame); + (void)WriteFrameDoubleDataset(file, root + "/nodal/triad", + {node_count, 3U, 3U}, triads.data(), triads.size(), + "X,Y,Z", "1,1,1", "global-cartesian", + "nodal-triad", frame); + const auto directors = FlattenDirectors(frame); + (void)WriteFrameDoubleDataset(file, root + "/nodal/director", {node_count, 3U}, + directors.data(), directors.size(), "D1,D2,D3", + "1,1,1", "global-cartesian", "nodal-director", + frame); + const auto internal = CopyVector(frame.internal_force); + (void)WriteFrameDoubleDataset( + file, root + "/nodal/generalized_internal_force", nodal_dimensions, + internal.data(), internal.size(), "FI1,FI2,FI3,MI_PHI1,MI_PHI2,MI_PHI3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "nodal", frame); + const auto external = CopyVector(frame.external_force); + (void)WriteFrameDoubleDataset( + file, root + "/nodal/generalized_external_force", nodal_dimensions, + external.data(), external.size(), "FE1,FE2,FE3,ME_PHI1,ME_PHI2,ME_PHI3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "nodal", frame); + const auto residual = CopyVector(frame.residual); + (void)WriteFrameDoubleDataset( + file, root + "/nodal/generalized_residual", nodal_dimensions, + residual.data(), residual.size(), "RU1,RU2,RU3,RPHI1,RPHI2,RPHI3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "nodal", frame); + const auto reaction = SpatialReaction(frame); + (void)WriteFrameDoubleDataset( + file, root + "/nodal/reaction", nodal_dimensions, reaction.data(), + reaction.size(), "RF1,RF2,RF3,RM1,RM2,RM3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "nodal", frame); + + std::vector local_frames; + std::vector strains; + std::vector resultants; + std::vector stresses; + local_frames.reserve(frame.shell_results.size() * 9U); + strains.reserve(frame.shell_results.size() * 8U); + resultants.reserve(frame.shell_results.size() * 8U); + stresses.reserve(frame.shell_results.size() * kSectionPositions * 3U); + for (const auto& row : frame.shell_results) { + for (const auto& axis : row.local_frame) { + local_frames.insert(local_frames.end(), axis.begin(), axis.end()); + } + strains.insert(strains.end(), row.generalized_strain.begin(), + row.generalized_strain.end()); + resultants.insert(resultants.end(), row.section_resultant.begin(), + row.section_resultant.end()); + for (const auto& stress : row.stress) { + stresses.insert(stresses.end(), stress.components.begin(), + stress.components.end()); + } + } + const hsize_t element_count = static_cast(domain.ShellElements().Size()); + const std::string shell_root = root + "/element/shell"; + auto local_frame = WriteFrameDoubleDataset( + file, shell_root + "/local_frame", {element_count, 4U, 3U, 3U}, + local_frames.data(), local_frames.size(), "X,Y,Z", "1,1,1", + "global-cartesian", "shell-local-frame", frame); + hdf5_internal::WriteStringAttribute(local_frame.Get(), "axis_names", "E1,E2,E3"); + hdf5_internal::WriteStringAttribute(local_frame.Get(), "source_element_type_dataset", + "/model/elements.source_element_type"); + hdf5_internal::WriteStringAttribute(local_frame.Get(), "internal_formulation", + "FESA-MITC4-total-lagrangian"); + hdf5_internal::WriteStringAttribute(local_frame.Get(), "midsurface_location_dataset", + "/model/shell/midsurface_locations"); + auto strain = WriteFrameDoubleDataset( + file, shell_root + "/generalized_strain", {element_count, 4U, 8U}, + strains.data(), strains.size(), "E11,E22,G12,K11,K22,K12,G13,G23", + "1,1,1,1/length,1/length,1/length,1,1", "shell-local", "midsurface", + frame); + WriteShellIdentity(strain.Get(), false, shell_root + "/local_frame"); + auto resultant = WriteFrameDoubleDataset( + file, shell_root + "/section_resultant", {element_count, 4U, 8U}, + resultants.data(), resultants.size(), "N11,N22,N12,M11,M22,M12,Q13,Q23", + "force/length,force/length,force/length,force,force,force,force/length," + "force/length", + "shell-local", "midsurface", frame); + WriteShellIdentity(resultant.Get(), false, shell_root + "/local_frame"); + auto stress = WriteFrameDoubleDataset( + file, shell_root + "/stress", {element_count, 4U, 3U, 3U}, stresses.data(), + stresses.size(), "S11,S22,S12", + "force/length^2,force/length^2,force/length^2", "shell-local", + "section-position", frame); + WriteShellIdentity(stress.Get(), true, shell_root + "/local_frame"); + (void)WriteFrameDoubleDataset(file, root + "/global/energy", {1U}, + &frame.physical_strain_energy, 1U, + "PHYSICAL_STRAIN_ENERGY", "force*length", + "global", "global", frame); + (void)WriteFrameDoubleDataset( + file, root + "/global/equilibrium", {6U}, frame.global_equilibrium.data(), + frame.global_equilibrium.size(), + "FORCE_1,FORCE_2,FORCE_3,MOMENT_1,MOMENT_2,MOMENT_3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "global-origin", frame); +} + +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(hdf5_internal::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 = hdf5_internal::MakeUtf8StringType(); + Hdf5Handle file_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(DiagnosticWriteRow)), + "Unable to create nonlinear diagnostic file type."), + H5Tclose}; + Hdf5Handle memory_type{hdf5_internal::RequireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(DiagnosticWriteRow)), + "Unable to create nonlinear diagnostic memory type."), + H5Tclose}; + const auto insert = [&](const hid_t type, const hid_t integer_type) { + hdf5_internal::RequireHdf5( + H5Tinsert(type, "severity", HOFFSET(DiagnosticWriteRow, severity), + string_type.Get()), + "Unable to define nonlinear diagnostic severity."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "code", HOFFSET(DiagnosticWriteRow, code), + string_type.Get()), + "Unable to define nonlinear diagnostic code."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "file", HOFFSET(DiagnosticWriteRow, file), + string_type.Get()), + "Unable to define nonlinear diagnostic file."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "line", HOFFSET(DiagnosticWriteRow, line), integer_type), + "Unable to define nonlinear diagnostic line."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "keyword", HOFFSET(DiagnosticWriteRow, keyword), + string_type.Get()), + "Unable to define nonlinear diagnostic keyword."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "entity_identity", + HOFFSET(DiagnosticWriteRow, entity_identity), string_type.Get()), + "Unable to define nonlinear diagnostic entity identity."); + hdf5_internal::RequireHdf5( + H5Tinsert(type, "message", HOFFSET(DiagnosticWriteRow, message), + string_type.Get()), + "Unable to define nonlinear diagnostic message."); + }; + insert(file_type.Get(), H5T_STD_U64LE); + insert(memory_type.Get(), H5T_NATIVE_UINT64); + auto dataset = hdf5_internal::WriteCompoundDataset( + file, "/diagnostics", rows.size(), file_type.Get(), memory_type.Get(), + rows.data()); + hdf5_internal::WriteUint64Attribute(dataset.Get(), "schema_version", 1U); +} + +Status ValidateDiagnostics(const std::vector& diagnostics) { + for (const auto& diagnostic : diagnostics) { + if (!hdf5_internal::IsValidUtf8(diagnostic.code) || + !hdf5_internal::IsValidUtf8(diagnostic.keyword) || + !hdf5_internal::IsValidUtf8(diagnostic.entity_identity) || + !hdf5_internal::IsValidUtf8(diagnostic.message)) { + return Failure("invalid-result-diagnostic", + "Nonlinear diagnostic text must be valid UTF-8."); + } + } + return Status::Ok(); +} + +} // namespace + +Status NonlinearHdf5ResultsWriter::Write( + const std::filesystem::path& output_path, const Domain& domain, + const NonlinearAnalysisState& state, + const std::vector& diagnostics) { + if (output_path.empty() || output_path.filename().empty()) { + return Failure("invalid-output-path", "The HDF5 output path must name a file."); + } + const Status domain_status = ValidateDomain(domain); + if (!domain_status.IsOk()) { + return domain_status; + } + const Status state_status = ValidateState(domain, state); + if (!state_status.IsOk()) { + return state_status; + } + const Status diagnostic_status = ValidateDiagnostics(diagnostics); + if (!diagnostic_status.IsOk()) { + return diagnostic_status; + } + NonlinearModelData model_data; + const Status model_status = BuildModelData(domain, model_data); + if (!model_status.IsOk()) { + return model_status; + } + + try { + hdf5_internal::Hdf5ErrorSilencer silence_backend_errors; + const std::filesystem::path candidate_path = + hdf5_internal::MakeCandidatePath(output_path); + hdf5_internal::CandidateFileGuard cleanup{candidate_path}; + Hdf5Handle file{hdf5_internal::RequireHdf5Id( + H5Fcreate(candidate_path.string().c_str(), H5F_ACC_EXCL, + H5P_DEFAULT, H5P_DEFAULT), + "Unable to create nonlinear temporary HDF5 file."), + H5Fclose}; + (void)hdf5_internal::CreateGroup(file.Get(), "/model"); + (void)hdf5_internal::CreateGroup(file.Get(), "/model/shell"); + (void)hdf5_internal::CreateGroup(file.Get(), "/steps/Step-1"); + WriteMetadata(file.Get(), domain); + WriteNodes(file.Get(), domain); + WriteElements(file.Get(), domain); + WriteModelData(file.Get(), domain, model_data); + WriteAcceptedFrameTable(file.Get(), state.AcceptedFrames()); + for (const auto& frame : state.AcceptedFrames()) { + WriteFrame(file.Get(), domain, frame); + } + WriteDiagnostics(file.Get(), diagnostics); + hdf5_internal::RequireHdf5(H5Fflush(file.Get(), H5F_SCOPE_GLOBAL), + "Unable to flush nonlinear temporary HDF5 file."); + hdf5_internal::RequireHdf5( + file.CloseChecked(), "Unable to close nonlinear temporary HDF5 file."); + const Status self_check = hdf5_internal::NonlinearSelfCheck( + candidate_path, domain, state, diagnostics.size()); + const Status atomic = hdf5_internal::AtomicFinalizeValidatedCandidate( + candidate_path, output_path, self_check); + if (atomic.IsOk()) { + cleanup.Release(); + } + return atomic; + } catch (const Hdf5Failure& failure) { + return Failure("hdf5-write-failure", failure.what()); + } catch (const std::exception& failure) { + return Failure("hdf5-write-failure", failure.what()); + } +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 663481d..43fc174 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -36,6 +36,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/nonlinear_hdf5_results_writer_test.cpp unit/io/hdf5/hdf5_writer_components_test.cpp unit/loads/load_test.cpp unit/model/domain_test.cpp diff --git a/tests/unit/io/hdf5/nonlinear_hdf5_results_writer_test.cpp b/tests/unit/io/hdf5/nonlinear_hdf5_results_writer_test.cpp new file mode 100644 index 0000000..b2cfb77 --- /dev/null +++ b/tests/unit/io/hdf5/nonlinear_hdf5_results_writer_test.cpp @@ -0,0 +1,591 @@ +#include "fesa/io/hdf5/nonlinear_hdf5_results_writer.h" + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/analysis/nonlinear_analysis_state.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/model/domain.h" +#include "fesa/results/nonlinear_results_writer.h" +#include "io/hdf5/hdf5_raii.h" + +namespace { + +using fesa::hdf5_internal::Hdf5Handle; + +constexpr std::array, 3> kIdentityTriad{ + std::array{1.0, 0.0, 0.0}, + std::array{0.0, 1.0, 0.0}, + std::array{0.0, 0.0, 1.0}}; + +class TempDirectory { + public: + explicit TempDirectory(const std::string& leaf) + : path_{std::filesystem::temp_directory_path() / + ("fesa-nonlinear-hdf5-" + leaf)} { + std::error_code error; + std::filesystem::remove_all(path_, error); + if (!std::filesystem::create_directory(path_)) { + throw std::runtime_error{"Unable to create nonlinear HDF5 test directory."}; + } + } + + ~TempDirectory() { + std::error_code error; + std::filesystem::remove_all(path_, error); + } + + const std::filesystem::path& Path() const noexcept { return path_; } + + private: + std::filesystem::path path_; +}; + +class WinHandle { + public: + explicit WinHandle(const HANDLE handle) : handle_{handle} {} + WinHandle(const WinHandle&) = delete; + WinHandle& operator=(const WinHandle&) = delete; + ~WinHandle() { + if (handle_ != INVALID_HANDLE_VALUE) { + (void)CloseHandle(handle_); + } + } + + HANDLE Get() const noexcept { return handle_; } + + private: + HANDLE handle_{INVALID_HANDLE_VALUE}; +}; + +struct Fixture { + std::unique_ptr domain; + std::unique_ptr dofs; + fesa::NonlinearAnalysisState state; +}; + +struct AcceptedFrameRow { + std::uint64_t frame_index; + std::uint64_t increment_index; + double load_factor; + std::uint64_t iteration_count; + std::uint8_t is_factor_one; +}; + +struct DiagnosticRow { + char* severity; + char* code; + char* file; + std::uint64_t line; + char* keyword; + char* entity_identity; + char* message; +}; + +fesa::ModelDefinition MakeDefinition(const std::filesystem::path& source) { + fesa::ModelDefinition definition{}; + definition.source_path = source; + definition.source_content_identity = "nonlinear-hdf5-fixture"; + definition.nodes = { + {{"Shell-1", 10, "10"}, {0.0, 0.0, 0.0}, {source, 1U}}, + {{"Shell-1", 20, "20"}, {1.0, 0.0, 0.0}, {source, 2U}}, + {{"Shell-1", 30, "30"}, {1.0, 1.0, 0.0}, {source, 3U}}, + {{"Shell-1", 40, "40"}, {0.0, 1.0, 0.0}, {source, 4U}}}; + definition.materials = {{"Steel", 210.0e9, 0.3, {source, 5U}}}; + definition.shell_sections = {{"Shell", 0.1, 0U, {source, 6U}}}; + definition.shell_elements = {{{"Shell-1", 1, "1"}, + fesa::ShellSourceElementType::kS4, + {0U, 1U, 2U, 3U}, + 0U, + 0U, + {source, 7U}}}; + for (std::size_t node = 0U; node < definition.nodes.size(); ++node) { + definition.shell_node_initial_frames.push_back( + {static_cast(node), + {0.0, 0.0, 1.0}, + {1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}}); + } + definition.steps = {{"Step-1", + {}, + {}, + fesa::StaticProcedureKind::kNonlinearFullNewton, + true, + 0.5, + 1.0, + 0.001, + 1.0, + {source, 8U}}}; + return definition; +} + +fesa::NonlinearStateBaselineInput MakeBaseline(const fesa::Domain& domain) { + fesa::NonlinearStateBaselineInput baseline{}; + for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { + baseline.node_baselines.push_back( + {static_cast(node), domain.Nodes()[node].source_id, + kIdentityTriad, {0.0, 0.0, 1.0}}); + } + return baseline; +} + +fesa::ShellResultRow MakeShellRow(const fesa::ShellMidsurfaceLocation location, + const std::array coordinates, + const double seed) { + return {0U, + location, + coordinates, + kIdentityTriad, + {seed, seed + 1.0, seed + 2.0, seed + 3.0, seed + 4.0, + seed + 5.0, seed + 6.0, seed + 7.0}, + {seed + 10.0, seed + 11.0, seed + 12.0, seed + 13.0, + seed + 14.0, seed + 15.0, seed + 16.0, seed + 17.0}, + {{{fesa::ShellSectionPosition::kBottom, + -1.0, + {seed + 20.0, seed + 21.0, seed + 22.0}}, + {fesa::ShellSectionPosition::kMiddle, + 0.0, + {seed + 23.0, seed + 24.0, seed + 25.0}}, + {fesa::ShellSectionPosition::kTop, + 1.0, + {seed + 26.0, seed + 27.0, seed + 28.0}}}}}; +} + +std::vector MakeShellRows(const double seed) { + const double gauss = 1.0 / std::sqrt(3.0); + return { + MakeShellRow(fesa::ShellMidsurfaceLocation::kGp1, {-gauss, -gauss}, + seed), + MakeShellRow(fesa::ShellMidsurfaceLocation::kGp2, {gauss, -gauss}, + seed + 30.0), + MakeShellRow(fesa::ShellMidsurfaceLocation::kGp3, {gauss, gauss}, + seed + 60.0), + MakeShellRow(fesa::ShellMidsurfaceLocation::kGp4, {-gauss, gauss}, + seed + 90.0)}; +} + +fesa::AcceptedFrameState MakeFrame(const fesa::NonlinearTrialState& committed, + const std::uint32_t frame_index, + const std::uint32_t increment_index, + const double load_factor) { + fesa::AcceptedFrameState frame{ + frame_index, + increment_index, + load_factor, + 2U, + load_factor == 1.0, + committed.nodal_internal_order, + committed.nodal_source_ids, + committed.generalized_coordinates, + committed.displacement, + committed.rotation_vector, + committed.nodal_triad, + committed.nodal_director, + committed.internal_force, + committed.external_force, + committed.residual, + committed.reaction, + MakeShellRows(100.0 * load_factor), + 20.0 * load_factor, + {load_factor, 2.0 * load_factor, 3.0 * load_factor, + 4.0 * load_factor, 5.0 * load_factor, 6.0 * load_factor}}; + for (std::size_t node = 0U; node < frame.nodal_internal_order.size(); ++node) { + const std::size_t offset = 6U * node; + frame.displacement[offset] = load_factor * static_cast(node + 1U); + frame.displacement[offset + 1U] = -frame.displacement[offset]; + frame.displacement[offset + 2U] = 0.5 * frame.displacement[offset]; + frame.rotation_vector[offset + 5U] = 0.1 * load_factor; + for (std::size_t component = 0U; component < 6U; ++component) { + frame.generalized_coordinates[offset + component] = + component < 3U ? frame.displacement[offset + component] + : frame.rotation_vector[offset + component]; + frame.internal_force[offset + component] = + 10.0 * static_cast(component + 1U); + frame.external_force[offset + component] = 0.0; + frame.residual[offset + component] = frame.internal_force[offset + component]; + frame.reaction[offset + component] = frame.residual[offset + component]; + } + } + return frame; +} + +Fixture MakeFixture(const std::filesystem::path& source) { + auto domain_result = fesa::Domain::Create(MakeDefinition(source)); + if (!domain_result.HasValue()) { + throw std::runtime_error{"Unable to construct nonlinear HDF5 Domain."}; + } + 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{"Unable to construct nonlinear HDF5 model."}; + } + auto dofs_result = fesa::DofManager::Create(model_result.Value()); + if (!dofs_result.HasValue()) { + throw std::runtime_error{"Unable to construct nonlinear HDF5 DOFs."}; + } + auto dofs = std::make_unique(std::move(dofs_result.Value())); + auto state_result = + fesa::NonlinearAnalysisState::Create(*domain, *dofs, MakeBaseline(*domain)); + if (!state_result.HasValue()) { + throw std::runtime_error{"Unable to construct nonlinear HDF5 state."}; + } + fesa::NonlinearAnalysisState state = std::move(state_result.Value()); + if (!state.CommitAcceptedFrame( + MakeFrame(state.CommittedState(), 0U, 1U, 0.5)) + .IsOk() || + !state.CommitAcceptedFrame( + MakeFrame(state.CommittedState(), 1U, 2U, 1.0)) + .IsOk()) { + throw std::runtime_error{"Unable to commit nonlinear HDF5 test frames."}; + } + return {std::move(domain), std::move(dofs), std::move(state)}; +} + +Hdf5Handle OpenFile(const std::filesystem::path& path) { + return {H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), H5Fclose}; +} + +Hdf5Handle OpenDataset(const hid_t file, const std::string& path) { + return {H5Dopen2(file, path.c_str(), H5P_DEFAULT), H5Dclose}; +} + +std::vector DatasetDimensions(const hid_t file, + const std::string& path) { + const auto dataset = OpenDataset(file, path); + Hdf5Handle space{H5Dget_space(dataset.Get()), H5Sclose}; + const int rank = H5Sget_simple_extent_ndims(space.Get()); + std::vector dimensions(static_cast(rank)); + if (rank > 0) { + EXPECT_GE(H5Sget_simple_extent_dims(space.Get(), dimensions.data(), nullptr), + 0); + } + return dimensions; +} + +std::string ReadStringAttribute(const hid_t object, const char* name) { + Hdf5Handle attribute{H5Aopen(object, name, H5P_DEFAULT), H5Aclose}; + Hdf5Handle type{H5Aget_type(attribute.Get()), H5Tclose}; + char* raw = nullptr; + EXPECT_GE(H5Aread(attribute.Get(), type.Get(), &raw), 0); + const std::string value = raw == nullptr ? std::string{} : std::string{raw}; + if (raw != nullptr) { + EXPECT_GE(H5free_memory(raw), 0); + } + return value; +} + +std::uint64_t ReadUint64Attribute(const hid_t object, const char* name) { + Hdf5Handle attribute{H5Aopen(object, name, H5P_DEFAULT), H5Aclose}; + std::uint64_t value = 0U; + EXPECT_GE(H5Aread(attribute.Get(), H5T_NATIVE_UINT64, &value), 0); + return value; +} + +double ReadDoubleAttribute(const hid_t object, const char* name) { + Hdf5Handle attribute{H5Aopen(object, name, H5P_DEFAULT), H5Aclose}; + double value = 0.0; + EXPECT_GE(H5Aread(attribute.Get(), H5T_NATIVE_DOUBLE, &value), 0); + return value; +} + +std::vector ReadDoubleDataset(const hid_t file, + const std::string& path) { + const auto dimensions = DatasetDimensions(file, path); + std::size_t count = 1U; + for (const hsize_t dimension : dimensions) { + count *= static_cast(dimension); + } + std::vector values(count); + const auto dataset = OpenDataset(file, path); + if (!values.empty()) { + EXPECT_GE(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), + 0); + } + return values; +} + +Hdf5Handle AcceptedFrameMemoryType() { + Hdf5Handle type{H5Tcreate(H5T_COMPOUND, sizeof(AcceptedFrameRow)), H5Tclose}; + EXPECT_GE(H5Tinsert(type.Get(), "frame_index", + HOFFSET(AcceptedFrameRow, frame_index), + H5T_NATIVE_UINT64), + 0); + EXPECT_GE(H5Tinsert(type.Get(), "increment_index", + HOFFSET(AcceptedFrameRow, increment_index), + H5T_NATIVE_UINT64), + 0); + EXPECT_GE(H5Tinsert(type.Get(), "load_factor", + HOFFSET(AcceptedFrameRow, load_factor), H5T_NATIVE_DOUBLE), + 0); + EXPECT_GE(H5Tinsert(type.Get(), "iteration_count", + HOFFSET(AcceptedFrameRow, iteration_count), + H5T_NATIVE_UINT64), + 0); + EXPECT_GE(H5Tinsert(type.Get(), "is_factor_one", + HOFFSET(AcceptedFrameRow, is_factor_one), + H5T_NATIVE_UINT8), + 0); + return type; +} + +std::vector ReadAcceptedFrames(const hid_t file) { + const auto dataset = + OpenDataset(file, "/steps/Step-1/accepted_frames"); + const auto dimensions = DatasetDimensions(file, "/steps/Step-1/accepted_frames"); + std::vector rows(static_cast(dimensions[0U])); + auto type = AcceptedFrameMemoryType(); + EXPECT_GE(H5Dread(dataset.Get(), type.Get(), H5S_ALL, H5S_ALL, H5P_DEFAULT, + rows.data()), + 0); + return rows; +} + +Hdf5Handle DiagnosticMemoryType() { + Hdf5Handle string_type{H5Tcopy(H5T_C_S1), H5Tclose}; + EXPECT_GE(H5Tset_size(string_type.Get(), H5T_VARIABLE), 0); + EXPECT_GE(H5Tset_cset(string_type.Get(), H5T_CSET_UTF8), 0); + Hdf5Handle type{H5Tcreate(H5T_COMPOUND, sizeof(DiagnosticRow)), H5Tclose}; + EXPECT_GE(H5Tinsert(type.Get(), "severity", HOFFSET(DiagnosticRow, severity), + string_type.Get()), + 0); + EXPECT_GE(H5Tinsert(type.Get(), "code", HOFFSET(DiagnosticRow, code), + string_type.Get()), + 0); + EXPECT_GE(H5Tinsert(type.Get(), "file", HOFFSET(DiagnosticRow, file), + string_type.Get()), + 0); + EXPECT_GE(H5Tinsert(type.Get(), "line", HOFFSET(DiagnosticRow, line), + H5T_NATIVE_UINT64), + 0); + EXPECT_GE(H5Tinsert(type.Get(), "keyword", HOFFSET(DiagnosticRow, keyword), + string_type.Get()), + 0); + EXPECT_GE(H5Tinsert(type.Get(), "entity_identity", + HOFFSET(DiagnosticRow, entity_identity), string_type.Get()), + 0); + EXPECT_GE(H5Tinsert(type.Get(), "message", HOFFSET(DiagnosticRow, message), + string_type.Get()), + 0); + return type; +} + +std::vector ReadDiagnostics(const hid_t file) { + const auto dataset = OpenDataset(file, "/diagnostics"); + const auto dimensions = DatasetDimensions(file, "/diagnostics"); + std::vector rows(static_cast(dimensions[0U])); + auto type = DiagnosticMemoryType(); + EXPECT_GE(H5Dread(dataset.Get(), type.Get(), H5S_ALL, H5S_ALL, H5P_DEFAULT, + rows.data()), + 0); + return rows; +} + +void ReclaimDiagnostics(const std::vector& rows) { + for (const auto& row : rows) { + for (char* value : {row.severity, row.code, row.file, row.keyword, + row.entity_identity, row.message}) { + if (value != nullptr) { + EXPECT_GE(H5free_memory(value), 0); + } + } + } +} + +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 nonlinear HDF5 test bytes."}; + } +} + +} // namespace + +// NSNR-HDF5-001 through NSNR-HDF5-006, NSNR-HDF5-008 through NSNR-HDF5-010. +TEST(NonlinearHdf5ResultsWriter, + WritesSchemaV1ModelAndCommittedAcceptedFramesOnly) { + TempDirectory directory{"schema"}; + auto fixture = MakeFixture(directory.Path() / "model.inp"); + const auto output = directory.Path() / "model_results.h5"; + + fesa::NonlinearHdf5ResultsWriter writer; + const auto status = writer.Write(output, *fixture.domain, fixture.state, {}); + ASSERT_TRUE(status.IsOk()) << status.Diagnostics().at(0U).message; + + const auto file = OpenFile(output); + ASSERT_GE(file.Get(), 0); + Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), H5Gclose}; + ASSERT_GE(metadata.Get(), 0); + EXPECT_EQ(ReadUint64Attribute(metadata.Get(), "schema_version"), 1U); + EXPECT_EQ(ReadStringAttribute(metadata.Get(), "feature_id"), + "nonlinear-static-newton-raphson"); + EXPECT_EQ(ReadStringAttribute(metadata.Get(), "analysis_procedure"), + "nonlinear-static-full-newton"); + EXPECT_EQ(ReadStringAttribute(metadata.Get(), "internal_formulation"), + "FESA-MITC4-total-lagrangian"); + EXPECT_EQ(ReadStringAttribute(metadata.Get(), "integration_rule"), + "2x2x2-gauss; mitc4-edge-midpoint-shear"); + EXPECT_EQ(ReadUint64Attribute(metadata.Get(), "factor_zero_is_frame"), 0U); + EXPECT_DOUBLE_EQ(ReadDoubleAttribute(metadata.Get(), "factor_zero_load_factor"), + 0.0); + + for (const char* path : + {"/model/nodes", "/model/elements", "/model/shell/materials", + "/model/shell/sections", "/model/shell/nodal_director", + "/model/shell/nodal_frame", "/model/nodal_constraint_mask", + "/model/prescribed_displacement_start", + "/model/prescribed_displacement_target", + "/model/shell/midsurface_locations", "/model/shell/section_positions"}) { + EXPECT_GT(H5Lexists(file.Get(), path, H5P_DEFAULT), 0) << path; + } + EXPECT_EQ(DatasetDimensions(file.Get(), "/steps/Step-1/accepted_frames"), + std::vector({2U})); + const auto frames = ReadAcceptedFrames(file.Get()); + ASSERT_EQ(frames.size(), 2U); + EXPECT_EQ(frames[0U].frame_index, 0U); + EXPECT_EQ(frames[0U].increment_index, 1U); + EXPECT_DOUBLE_EQ(frames[0U].load_factor, 0.5); + EXPECT_EQ(frames[0U].is_factor_one, 0U); + EXPECT_EQ(frames[1U].frame_index, 1U); + EXPECT_EQ(frames[1U].increment_index, 2U); + EXPECT_DOUBLE_EQ(frames[1U].load_factor, 1.0); + EXPECT_EQ(frames[1U].is_factor_one, 1U); + EXPECT_EQ(H5Lexists(file.Get(), "/steps/Step-1/frames/2", H5P_DEFAULT), 0); + EXPECT_EQ(H5Lexists(file.Get(), "/steps/Step-1/trial_history", H5P_DEFAULT), + 0); + + const std::string root = "/steps/Step-1/frames/1"; + for (const std::string& suffix : + {"/nodal/displacement", "/nodal/triad", "/nodal/director", + "/nodal/generalized_internal_force", "/nodal/generalized_external_force", + "/nodal/generalized_residual", "/nodal/reaction", + "/element/shell/local_frame", "/element/shell/generalized_strain", + "/element/shell/section_resultant", "/element/shell/stress", + "/global/energy", "/global/equilibrium"}) { + EXPECT_GT(H5Lexists(file.Get(), (root + suffix).c_str(), H5P_DEFAULT), 0) + << suffix; + } + EXPECT_EQ(DatasetDimensions(file.Get(), root + "/nodal/displacement"), + std::vector({4U, 6U})); + EXPECT_EQ(DatasetDimensions(file.Get(), root + "/nodal/triad"), + std::vector({4U, 3U, 3U})); + EXPECT_EQ(DatasetDimensions(file.Get(), root + "/element/shell/local_frame"), + std::vector({1U, 4U, 3U, 3U})); + EXPECT_EQ(DatasetDimensions(file.Get(), root + "/element/shell/stress"), + std::vector({1U, 4U, 3U, 3U})); + const auto displacement = ReadDoubleDataset(file.Get(), root + "/nodal/displacement"); + ASSERT_EQ(displacement.size(), 24U); + EXPECT_DOUBLE_EQ(displacement[0U], 1.0); + EXPECT_DOUBLE_EQ(displacement[1U], -1.0); + EXPECT_DOUBLE_EQ(displacement[2U], 0.5); + EXPECT_DOUBLE_EQ(displacement[5U], 0.1); +} + +// NSNR-HDF5-007. +TEST(NonlinearHdf5ResultsWriter, ProjectsRotationalReactionToSpatialMoment) { + TempDirectory directory{"reaction"}; + auto fixture = MakeFixture(directory.Path() / "reaction.inp"); + const auto output = directory.Path() / "reaction_results.h5"; + + fesa::NonlinearHdf5ResultsWriter writer; + ASSERT_TRUE(writer.Write(output, *fixture.domain, fixture.state, {}).IsOk()); + + const auto file = OpenFile(output); + const auto reaction = ReadDoubleDataset( + file.Get(), "/steps/Step-1/frames/1/nodal/reaction"); + ASSERT_EQ(reaction.size(), 24U); + const double theta = 0.1; + const double quadratic = + (1.0 - 0.5 * theta / std::tan(0.5 * theta)) / (theta * theta); + EXPECT_NEAR(reaction[3U], + (1.0 - quadratic * theta * theta) * 40.0 - + 0.5 * theta * 50.0, + 1.0e-12); + EXPECT_NEAR(reaction[4U], + 0.5 * theta * 40.0 + + (1.0 - quadratic * theta * theta) * 50.0, + 1.0e-12); + EXPECT_NEAR(reaction[5U], 60.0, 1.0e-12); + EXPECT_GT(quadratic, 0.0); +} + +// NSNR-HDF5-012 and NSNR-HDF5-013. +TEST(NonlinearHdf5ResultsWriter, StoresSortedDiagnosticsThroughBackendNeutralBoundary) { + static_assert(std::is_base_of_v); + TempDirectory directory{"diagnostics"}; + auto fixture = MakeFixture(directory.Path() / "diagnostics.inp"); + const std::vector diagnostics{ + {fesa::Severity::kWarning, + "later", + {fixture.domain->SourcePath(), 20U}, + "*OUTPUT", + "", + "Later warning."}, + {fesa::Severity::kWarning, + "earlier", + {fixture.domain->SourcePath(), 10U}, + "*PREPRINT", + "", + "Earlier warning."}}; + const auto output = directory.Path() / "diagnostics_results.h5"; + + fesa::NonlinearHdf5ResultsWriter writer; + ASSERT_TRUE(writer.Write(output, *fixture.domain, fixture.state, diagnostics) + .IsOk()); + const auto file = OpenFile(output); + const auto rows = ReadDiagnostics(file.Get()); + ASSERT_EQ(rows.size(), 2U); + EXPECT_STREQ(rows[0U].code, "earlier"); + EXPECT_EQ(rows[0U].line, 10U); + EXPECT_STREQ(rows[1U].code, "later"); + EXPECT_EQ(rows[1U].line, 20U); + ReclaimDiagnostics(rows); +} + +// NSNR-HDF5-011. +TEST(NonlinearHdf5ResultsWriter, CandidateFinalizationFailurePreservesPriorFinal) { + TempDirectory directory{"atomic"}; + auto fixture = MakeFixture(directory.Path() / "atomic.inp"); + const auto output = directory.Path() / "atomic_results.h5"; + const std::vector sentinel{'p', 'r', 'i', 'o', 'r'}; + WriteBytes(output, sentinel); + WinHandle lock{CreateFileW(output.c_str(), GENERIC_READ, + FILE_SHARE_READ | FILE_SHARE_WRITE, nullptr, + OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, nullptr)}; + ASSERT_NE(lock.Get(), INVALID_HANDLE_VALUE); + + fesa::NonlinearHdf5ResultsWriter writer; + const auto status = writer.Write(output, *fixture.domain, fixture.state, {}); + + ASSERT_FALSE(status.IsOk()); + EXPECT_EQ(status.Category(), fesa::FailureCategory::kOutput); + ASSERT_EQ(status.Diagnostics().size(), 1U); + EXPECT_EQ(status.Diagnostics()[0U].code, "hdf5-finalization-failure"); + EXPECT_EQ(ReadBytes(output), sentinel); +}