diff --git a/docs/HDF5_SCHEMA.md b/docs/HDF5_SCHEMA.md index 90e2c6d..fbd9e0d 100644 --- a/docs/HDF5_SCHEMA.md +++ b/docs/HDF5_SCHEMA.md @@ -1,90 +1,164 @@ -# FESA HDF5 Schema 1.0.0 +# FESA HDF5 Schema 2.0.0 -## 1. 범위 +## 1. Scope and version compatibility -Schema `1.0.0`은 `results-and-pipeline` Phase의 최소 수직 슬라이스를 정의한다. -파일은 활성 `Domain`의 절점, Beam 연결성, 적용된 전단면적과 그 출처, 그리고 전역 -좌표계 절점 변위·회전 및 반력·반력모멘트를 저장한다. 단위 변환은 수행하지 않는다. +Schema `2.0.0` is the self-contained Phase 1 result contract. One file contains +the active normalized model, its single linear-static analysis definition and +solver settings, and every nodal and Beam result needed without the source +`.inp` file. FESA performs no unit conversion. -이 버전에는 재료 전체 속성, 집합, 하중·경계조건, solver 설정, 요소 결과, history, -reference CSV 및 진단 dataset을 저장하지 않는다. 이후 같은 major version에서 -dataset을 추가할 수 있지만 아래 required object의 의미, 형상 또는 datatype을 -변경해서는 안 된다. +Schema `1.0.0` was the earlier minimal vertical slice. Version `2.0.0` changes +the required model and result objects, so it is a new major version rather than +an in-place change to `1.0.0`. The current writer and reader accept exactly +`2.0.0`; every other version fails with `hdf5.unsupported_schema`. A future +reader may explicitly add support for compatible minor versions, but must not +infer compatibility from a version prefix. -## 2. 공통 규칙 +## 2. Common rules -- root attribute `schema_version`은 UTF-8 문자열 `1.0.0`이다. -- 정수 dataset은 명시한 little-endian 고정폭 타입을 사용한다. -- 실수 dataset은 IEEE 754 little-endian 64-bit 타입을 사용한다. -- 문자열 dataset과 attribute는 UTF-8 variable-length string을 사용한다. -- `dense_index`는 해당 dataset 행의 0-based index이며 연속적이다. -- `internal_id`와 결과의 `node_ids`는 FESA semantic model의 nonnegative ID다. -- flat/orphan mesh의 `part_name`과 `instance_name`은 빈 문자열이다. -- 결과의 6개 component 순서는 - `(Ux, Uy, Uz, Rx, Ry, Rz)` 및 `(RFx, RFy, RFz, RMx, RMy, RMz)`다. -- Step과 frame group 이름은 각각 0부터 연속된 decimal index다. 원래 Step 이름은 - Step group의 `name` attribute에 저장한다. +- Root attributes are variable-length UTF-8 strings: + `schema_version="2.0.0"`, `fesa_version`, and + `unit_policy="consistent_input_units_no_conversion"`. +- Integer datasets use the stated little-endian fixed-width type. Floating + datasets use IEEE 754 little-endian `float64`. Strings are variable-length + UTF-8. +- `dense_index` is a contiguous 0-based row index. Semantic `internal_id` + values are nonnegative and are not assumed to be dense or ordered. +- Flat/orphan mesh `part_name` and `instance_name` values are empty strings. +- Ragged arrays use an offset dataset of length `row_count + 1`. Offsets start + at zero, are nondecreasing, and the final offset equals the flattened row + count. Input order is preserved. +- Numeric field datasets carry UTF-8 `coordinate_system` and `components` + attributes where listed. `components` is a comma-separated ordered list. +- Step and frame group names are contiguous decimal indices beginning at zero. + Phase 1 requires exactly one analysis step and one result step with the same + name. ## 3. Required objects -```text -/ -├── @schema_version UTF-8 = "1.0.0" -├── model -│ ├── nodes -│ │ ├── dense_index uint64 [node_count] -│ │ ├── internal_id int64 [node_count] -│ │ ├── part_name UTF-8 [node_count] -│ │ ├── instance_name UTF-8 [node_count] -│ │ ├── local_label int64 [node_count] -│ │ └── coordinates float64[node_count, 3] -│ ├── elements -│ │ ├── dense_index uint64 [element_count] -│ │ ├── internal_id int64 [element_count] -│ │ ├── connectivity uint64 [element_count, 2] -│ │ └── section_id int64 [element_count] -│ └── sections -│ ├── internal_id int64 [section_count] -│ ├── shear_area_y float64[section_count] -│ ├── shear_area_z float64[section_count] -│ └── shear_source uint8 [section_count] -└── results - └── steps - └── - ├── @name UTF-8 - └── frames - └── - ├── @step_time float64 - └── nodal - ├── node_ids int64 [result_node_count] - ├── displacement float64[result_node_count, 6] - └── reaction float64[result_node_count, 6] -``` +### 3.1 Model -`model/elements/connectivity`는 `model/nodes/dense_index`를 참조한다. 따라서 -`internal_id`가 연속적이거나 Domain 저장 순서와 같다고 가정하지 않는다. -`section_id`는 `model/sections/internal_id`를 참조한다. +Let `N`, `E`, `M`, `S`, `NS`, and `ES` be the node, Beam element, material, +section, node-set, and element-set counts. Let `P` be the total number of +section recovery points, `NM` the total node-set membership count, and `EM` the +total element-set membership count. -`shear_source` 값은 다음과 같다. +| Path | Type | Rank and shape | Attributes / meaning | +|---|---|---|---| +| `/model/nodes/dense_index` | `uint64` | 1, `[N]` | contiguous row index | +| `/model/nodes/internal_id` | `int64` | 1, `[N]` | `NodeId` | +| `/model/nodes/part_name` | UTF-8 | 1, `[N]` | entity provenance | +| `/model/nodes/instance_name` | UTF-8 | 1, `[N]` | entity provenance | +| `/model/nodes/local_label` | `int64` | 1, `[N]` | external local label | +| `/model/nodes/coordinates` | `float64` | 2, `[N,3]` | `coordinate_system="global"`, `components="X,Y,Z"` | +| `/model/elements/dense_index` | `uint64` | 1, `[E]` | contiguous row index | +| `/model/elements/internal_id` | `int64` | 1, `[E]` | `ElementId` | +| `/model/elements/part_name` | UTF-8 | 1, `[E]` | entity provenance | +| `/model/elements/instance_name` | UTF-8 | 1, `[E]` | entity provenance | +| `/model/elements/local_label` | `int64` | 1, `[E]` | external local label | +| `/model/elements/connectivity` | `uint64` | 2, `[E,2]` | node `dense_index`, ordered end `-1,+1` | +| `/model/elements/material_id` | `int64` | 1, `[E]` | references material `internal_id` | +| `/model/elements/section_id` | `int64` | 1, `[E]` | references section `internal_id` | +| `/model/materials/internal_id` | `int64` | 1, `[M]` | `MaterialId` | +| `/model/materials/name` | UTF-8 | 1, `[M]` | material name | +| `/model/materials/young_modulus` | `float64` | 1, `[M]` | finite, positive | +| `/model/materials/poisson_ratio` | `float64` | 1, `[M]` | finite, `-1 < nu < 0.5` | +| `/model/sections/internal_id` | `int64` | 1, `[S]` | `SectionId` | +| `/model/sections/name` | UTF-8 | 1, `[S]` | section name | +| `/model/sections/area` | `float64` | 1, `[S]` | `A` | +| `/model/sections/moment_y` | `float64` | 1, `[S]` | `Iy` | +| `/model/sections/moment_z` | `float64` | 1, `[S]` | `Iz` | +| `/model/sections/torsion_constant` | `float64` | 1, `[S]` | `J` | +| `/model/sections/shear_area_y` | `float64` | 1, `[S]` | applied `Asy` | +| `/model/sections/shear_area_z` | `float64` | 1, `[S]` | applied `Asz` | +| `/model/sections/shear_source` | `uint8` | 1, `[S]` | `0=input`, `1=phase1_default` | +| `/model/sections/orientation` | `float64` | 2, `[S,3]` | `coordinate_system="global"`, `components="X,Y,Z"` | +| `/model/sections/recovery_point_offsets` | `uint64` | 1, `[S+1]` | offsets into `recovery_points` | +| `/model/sections/recovery_points` | `float64` | 2, `[P,2]` | `coordinate_system="element_local"`, `components="y,z"` | +| `/model/sets/node/names` | UTF-8 | 1, `[NS]` | exact node-set names | +| `/model/sets/node/member_offsets` | `uint64` | 1, `[NS+1]` | offsets into `members` | +| `/model/sets/node/members` | `int64` | 1, `[NM]` | `NodeId`, input set/member order | +| `/model/sets/element/names` | UTF-8 | 1, `[ES]` | exact element-set names | +| `/model/sets/element/member_offsets` | `uint64` | 1, `[ES+1]` | offsets into `members` | +| `/model/sets/element/members` | `int64` | 1, `[EM]` | `ElementId`, input set/member order | -| 값 | 의미 | +### 3.2 Analysis + +`/analysis/steps/0` has UTF-8 attribute `name`. Let `B` be the prescribed-DOF +count and `L` the nodal-load count. + +| Path | Type | Rank and shape | Attributes / meaning | +|---|---|---|---| +| `/analysis/steps/0/boundary_conditions/node_ids` | `int64` | 1, `[B]` | target `NodeId` | +| `/analysis/steps/0/boundary_conditions/dofs` | `uint8` | 1, `[B]` | Abaqus/FESA DOF number 1 through 6 | +| `/analysis/steps/0/boundary_conditions/values` | `float64` | 1, `[B]` | prescribed value | +| `/analysis/steps/0/nodal_loads/node_ids` | `int64` | 1, `[L]` | target `NodeId` | +| `/analysis/steps/0/nodal_loads/values` | `float64` | 2, `[L,6]` | `coordinate_system="global"`, `components="Fx,Fy,Fz,Mx,My,Mz"` | + +`/analysis/solver_settings` has the exact UTF-8 attributes used by the Phase 1 +pipeline: `backend="mkl_pardiso"`, +`matrix_storage="symmetric_upper_csr"`, +`matrix_type="symmetric_positive_definite"`, +`constraint_method="essential_dof_elimination"`, and +`assembly="deterministic_serial"`. No unused future settings are stored. + +### 3.3 Results + +`/results/steps/0` has UTF-8 attribute `name`; frame group `0` has scalar +`float64` attribute `step_time`. Let `RN`, `RE`, `RP`, and `D` be the nodal +result, Beam result, flattened Beam recovery-point, and diagnostic counts. + +| Path below `/results/steps/0/frames/0` | Type | Rank and shape | Attributes / meaning | +|---|---|---|---| +| `nodal/node_ids` | `int64` | 1, `[RN]` | `NodeId`; provenance is joined from `/model/nodes` | +| `nodal/displacement` | `float64` | 2, `[RN,6]` | `coordinate_system="global"`, `components="Ux,Uy,Uz,Rx,Ry,Rz"` | +| `nodal/reaction` | `float64` | 2, `[RN,6]` | `coordinate_system="global"`, `components="RFx,RFy,RFz,RMx,RMy,RMz"` | +| `element/beam/element_ids` | `int64` | 1, `[RE]` | `ElementId` | +| `element/beam/part_name` | UTF-8 | 1, `[RE]` | result provenance | +| `element/beam/instance_name` | UTF-8 | 1, `[RE]` | result provenance | +| `element/beam/local_label` | `int64` | 1, `[RE]` | result provenance | +| `element/beam/local_frame` | `float64` | 3, `[RE,3,3]` | `coordinate_system="global"`, `components="ex,ey,ez"`; last dimension is `X,Y,Z` | +| `element/beam/end_node_ids` | `int64` | 2, `[RE,2]` | ordered ends `-1,+1` | +| `element/beam/xi` | `float64` | 2, `[RE,2]` | `components="end_minus,end_plus"` | +| `element/beam/section_strain` | `float64` | 3, `[RE,2,6]` | `coordinate_system="element_local"`, `components="epsilon,gamma_y,gamma_z,kappa_x,kappa_y,kappa_z"` | +| `element/beam/section_force` | `float64` | 3, `[RE,2,6]` | `coordinate_system="element_local"`, `components="N,Vy,Vz,T,My,Mz"` | +| `element/beam/centroid_sigma_xx` | `float64` | 2, `[RE,2]` | `coordinate_system="element_local"`, `components="end_minus,end_plus"`, `quantity="sigma_xx"` | +| `element/beam/recovery_point_offsets` | `uint64` | 1, `[RE+1]` | offsets into recovery-point rows | +| `element/beam/recovery_point_sigma_xx` | `float64` | 2, `[RP,2]` | `coordinate_system="element_local"`, `components="end_minus,end_plus"`, `quantity="sigma_xx"`; point order comes from the referenced section | +| `diagnostics/stage` | `uint8` | 1, `[D]` | enum table below | +| `diagnostics/severity` | `uint8` | 1, `[D]` | `0=warning`, `1=error` | +| `diagnostics/code` | UTF-8 | 1, `[D]` | exact diagnostic code | +| `diagnostics/message` | UTF-8 | 1, `[D]` | exact diagnostic message | +| `diagnostics/has_source` | `uint8` | 1, `[D]` | `0=no source`, `1=source present` | +| `diagnostics/source_file` | UTF-8 | 1, `[D]` | empty when source is absent | +| `diagnostics/source_line` | `uint64` | 1, `[D]` | zero when source is absent | +| `diagnostics/source_column` | `uint64` | 1, `[D]` | zero when source is absent | + +Diagnostic stage encoding follows the declaration order: + +| Value | Stage | |---:|---| -| `0` | 입력에서 명시된 전단강성으로부터 구성한 값 (`input`) | -| `1` | Phase 1 기본값 `Asy=Asz=5A/6`, `SCF=0` (`phase1_default`) | +| 0 | `io` | +| 1 | `syntax` | +| 2 | `semantic` | +| 3 | `model` | +| 4 | `equation` | +| 5 | `solver` | +| 6 | `results` | +| 7 | `validation` | -## 4. Writer와 reader 계약 +## 4. Writer and reader contract -- writer는 쓰기 전에 `ResultDatabase` 유효성과 schema version을 검사한다. -- schema `1.0.0`이 표현하지 않는 non-empty frame diagnostics는 파일을 만들기 전에 - `hdf5.unsupported_result_diagnostics`로 거부한다. -- required object 생성·쓰기·flush·close 중 HDF5 오류가 발생하면 성공으로 반환하지 - 않고 `DiagnosticStage::results` 오류로 변환한다. -- reader는 schema version, required object, datatype, rank와 shape를 검사한다. -- reader는 model internal ID의 uniqueness, finite coordinates 및 finite positive - shear area를 검사한다. -- writer와 reader는 모든 nodal result ID가 `model/nodes/internal_id`에 존재하는지 - 검사하며, 없는 ID를 성공 결과로 반환하지 않는다. -- reader는 nodal result를 `ResultDatabase`로, model dataset을 HDF5 adapter 전용 - read-only inspection model로 반환한다. `Domain`과 `ResultDatabase`에는 HDF5 - 저장 계약을 추가하지 않는다. -- malformed 또는 지원하지 않는 파일은 부분 database를 반환하지 않는다. +- The writer validates `ResultDatabase`, exact schema version, model/result ID + and provenance joins, Beam connectivity, recovery-point counts, and the + single-step name before creating the file. +- Required-object creation, write, flush, and close failures become + `DiagnosticStage::results` errors. A failed write is never reported as + success. +- The reader validates the exact version, required datatypes, ranks, shapes, + offsets, finite values, uniqueness, references, field metadata, and result + contracts. It does not return a partial database or partial snapshots. +- The public reader returns adapter-owned, read-only metadata, model, and + analysis snapshots plus the semantic `ResultDatabase`. No HDF5 object or + handle escapes the adapter, and the snapshots contain enough information to + reconstruct the Phase 1 model and run definition without the source deck. diff --git a/include/fesa/io/hdf5/writer.hpp b/include/fesa/io/hdf5/writer.hpp index f7c0fd1..4139ebf 100644 --- a/include/fesa/io/hdf5/writer.hpp +++ b/include/fesa/io/hdf5/writer.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include @@ -26,27 +27,60 @@ struct Hdf5NodeSnapshot final { struct Hdf5ElementSnapshot final { std::uint64_t dense_index; ElementId id; + EntityOrigin origin; std::array connectivity; + MaterialId material; SectionId section; }; struct Hdf5SectionSnapshot final { SectionId id; + std::string name; + double area; + double iy; + double iz; + double torsion_j; double shear_area_y; double shear_area_z; ShearPropertySource shear_source; + Vec3 orientation; + std::vector> recovery_points; }; struct Hdf5ModelSnapshot final { std::vector nodes; std::vector elements; + std::vector materials; std::vector sections; + std::vector node_sets; + std::vector element_sets; +}; + +struct Hdf5MetadataSnapshot final { + std::string schema_version; + std::string fesa_version; + std::string unit_policy; +}; + +struct Hdf5SolverSettingsSnapshot final { + std::string backend; + std::string matrix_storage; + std::string matrix_type; + std::string constraint_method; + std::string assembly; +}; + +struct Hdf5AnalysisSnapshot final { + StepDefinition step; + Hdf5SolverSettingsSnapshot solver; }; struct Hdf5ReadResult final { std::optional database; std::vector diagnostics; std::optional model; + std::optional metadata; + std::optional analysis; }; [[nodiscard]] std::vector write_hdf5( diff --git a/src/fesa/analysis/linear_static_analysis.cpp b/src/fesa/analysis/linear_static_analysis.cpp index 6aa0e0a..9ef1632 100644 --- a/src/fesa/analysis/linear_static_analysis.cpp +++ b/src/fesa/analysis/linear_static_analysis.cpp @@ -203,7 +203,7 @@ AnalysisRunResult LinearStaticAnalysis::run(const Domain& domain) const { } ResultDatabase database{ - "1.0.0", + "2.0.0", {{ domain.step().name, {{ diff --git a/src/fesa/io/hdf5/writer.cpp b/src/fesa/io/hdf5/writer.cpp index 783ccee..3cc444a 100644 --- a/src/fesa/io/hdf5/writer.cpp +++ b/src/fesa/io/hdf5/writer.cpp @@ -1,5 +1,7 @@ #include +#include + #include #include @@ -21,7 +23,16 @@ namespace fesa { namespace { -constexpr std::string_view schema_version = "1.0.0"; +constexpr std::string_view schema_version = "2.0.0"; +constexpr std::string_view unit_policy = + "consistent_input_units_no_conversion"; +constexpr std::string_view solver_backend = "mkl_pardiso"; +constexpr std::string_view solver_matrix_storage = "symmetric_upper_csr"; +constexpr std::string_view solver_matrix_type = + "symmetric_positive_definite"; +constexpr std::string_view solver_constraint_method = + "essential_dof_elimination"; +constexpr std::string_view solver_assembly = "deterministic_serial"; class Hdf5Error final : public std::runtime_error { public: @@ -424,7 +435,10 @@ void write_numeric_2d( const hid_t memory_type, const std::size_t row_count, const std::size_t column_count, - const std::span values) { + const std::span values, + const std::string_view coordinate_system = {}, + const std::string_view components = {}, + const std::string_view quantity = {}) { if (row_count * column_count != values.size()) { fail( "hdf5.write_failed", @@ -465,6 +479,81 @@ void write_numeric_2d( "write dataset '" + owned_name + "'", "hdf5.write_failed"); } + if (!coordinate_system.empty()) { + write_string_attribute( + context, dataset.get(), "coordinate_system", coordinate_system); + } + if (!components.empty()) { + write_string_attribute( + context, dataset.get(), "components", components); + } + if (!quantity.empty()) { + write_string_attribute(context, dataset.get(), "quantity", quantity); + } +} + +template +void write_numeric_3d( + Hdf5Context& context, + const hid_t parent, + const std::string_view name, + const hid_t file_type, + const hid_t memory_type, + const std::size_t first_count, + const std::size_t second_count, + const std::size_t third_count, + const std::span values, + const std::string_view coordinate_system = {}, + const std::string_view components = {}) { + if (first_count * second_count * third_count != values.size()) { + fail( + "hdf5.write_failed", + "Internal HDF5 dataset shape does not match its values."); + } + const std::array dimensions{ + static_cast(first_count), + static_cast(second_count), + static_cast(third_count), + }; + auto space = create_space(context, dimensions, "hdf5.write_failed"); + const std::string owned_name{name}; + auto dataset = checked_handle( + context, + call_hdf5([&] { + return H5Dcreate2( + parent, + owned_name.c_str(), + file_type, + space.get(), + H5P_DEFAULT, + H5P_DEFAULT, + H5P_DEFAULT); + }), + &H5Dclose, + "dataset '" + owned_name + "'", + "hdf5.write_failed"); + if (!values.empty()) { + check_status( + call_hdf5([&] { + return H5Dwrite( + dataset.get(), + memory_type, + H5S_ALL, + H5S_ALL, + H5P_DEFAULT, + values.data()); + }), + "write dataset '" + owned_name + "'", + "hdf5.write_failed"); + } + if (!coordinate_system.empty()) { + write_string_attribute( + context, dataset.get(), "coordinate_system", coordinate_system); + } + if (!components.empty()) { + write_string_attribute( + context, dataset.get(), "components", components); + } } void write_strings_1d( @@ -582,16 +671,26 @@ void write_model( H5T_NATIVE_DOUBLE, domain.nodes().size(), 3U, - std::span{coordinates}); + std::span{coordinates}, + "global", + "X,Y,Z"); } std::vector element_dense; std::vector element_ids; + std::vector element_part_names; + std::vector element_instance_names; + std::vector element_local_labels; std::vector connectivity; + std::vector element_material_ids; std::vector element_section_ids; element_dense.reserve(domain.beam_elements().size()); element_ids.reserve(domain.beam_elements().size()); + element_part_names.reserve(domain.beam_elements().size()); + element_instance_names.reserve(domain.beam_elements().size()); + element_local_labels.reserve(domain.beam_elements().size()); connectivity.reserve(domain.beam_elements().size() * 2U); + element_material_ids.reserve(domain.beam_elements().size()); element_section_ids.reserve(domain.beam_elements().size()); for (std::size_t index = 0; index < domain.beam_elements().size(); ++index) { const BeamElement& element = domain.beam_elements()[index]; @@ -605,8 +704,12 @@ void write_model( } element_dense.push_back(static_cast(index)); element_ids.push_back(element.id.value()); + element_part_names.push_back(element.origin.part_name); + element_instance_names.push_back(element.origin.instance_name); + element_local_labels.push_back(element.origin.local_label); connectivity.push_back(first->second); connectivity.push_back(second->second); + element_material_ids.push_back(element.material.value()); element_section_ids.push_back(element.section.value()); } @@ -627,6 +730,23 @@ void write_model( H5T_STD_I64LE, H5T_NATIVE_INT64, std::span{element_ids}); + write_strings_1d( + context, + elements_group.get(), + "part_name", + element_part_names); + write_strings_1d( + context, + elements_group.get(), + "instance_name", + element_instance_names); + write_numeric_1d( + context, + elements_group.get(), + "local_label", + H5T_STD_I64LE, + H5T_NATIVE_INT64, + std::span{element_local_labels}); write_numeric_2d( context, elements_group.get(), @@ -636,6 +756,13 @@ void write_model( domain.beam_elements().size(), 2U, std::span{connectivity}); + write_numeric_1d( + context, + elements_group.get(), + "material_id", + H5T_STD_I64LE, + H5T_NATIVE_INT64, + std::span{element_material_ids}); write_numeric_1d( context, elements_group.get(), @@ -645,20 +772,90 @@ void write_model( std::span{element_section_ids}); } + { + auto materials_group = + create_group(context, model_group.get(), "materials"); + std::vector ids; + std::vector names; + std::vector young; + std::vector poisson; + ids.reserve(domain.materials().size()); + names.reserve(domain.materials().size()); + young.reserve(domain.materials().size()); + poisson.reserve(domain.materials().size()); + for (const IsotropicElastic& material : domain.materials()) { + ids.push_back(material.id.value()); + names.push_back(material.name); + young.push_back(material.young); + poisson.push_back(material.poisson); + } + write_numeric_1d( + context, + materials_group.get(), + "internal_id", + H5T_STD_I64LE, + H5T_NATIVE_INT64, + std::span{ids}); + write_strings_1d(context, materials_group.get(), "name", names); + write_numeric_1d( + context, + materials_group.get(), + "young_modulus", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + std::span{young}); + write_numeric_1d( + context, + materials_group.get(), + "poisson_ratio", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + std::span{poisson}); + } + std::vector section_ids; + std::vector section_names; + std::vector areas; + std::vector moment_y; + std::vector moment_z; + std::vector torsion_constants; std::vector shear_area_y; std::vector shear_area_z; std::vector shear_sources; + std::vector orientations; + std::vector recovery_offsets{0U}; + std::vector recovery_points; section_ids.reserve(domain.sections().size()); + section_names.reserve(domain.sections().size()); + areas.reserve(domain.sections().size()); + moment_y.reserve(domain.sections().size()); + moment_z.reserve(domain.sections().size()); + torsion_constants.reserve(domain.sections().size()); shear_area_y.reserve(domain.sections().size()); shear_area_z.reserve(domain.sections().size()); shear_sources.reserve(domain.sections().size()); + orientations.reserve(domain.sections().size() * 3U); + recovery_offsets.reserve(domain.sections().size() + 1U); for (const BeamSection& section : domain.sections()) { section_ids.push_back(section.id.value()); + section_names.push_back(section.name); + areas.push_back(section.area); + moment_y.push_back(section.iy); + moment_z.push_back(section.iz); + torsion_constants.push_back(section.torsion_j); shear_area_y.push_back(section.shear_area_y); shear_area_z.push_back(section.shear_area_z); shear_sources.push_back( section.shear_source == ShearPropertySource::input ? 0U : 1U); + orientations.push_back(section.orientation.x); + orientations.push_back(section.orientation.y); + orientations.push_back(section.orientation.z); + for (const auto& point : section.recovery_points) { + recovery_points.push_back(point[0]); + recovery_points.push_back(point[1]); + } + recovery_offsets.push_back( + static_cast(recovery_points.size() / 2U)); } { @@ -671,6 +868,36 @@ void write_model( H5T_STD_I64LE, H5T_NATIVE_INT64, std::span{section_ids}); + write_strings_1d( + context, sections_group.get(), "name", section_names); + write_numeric_1d( + context, + sections_group.get(), + "area", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + std::span{areas}); + write_numeric_1d( + context, + sections_group.get(), + "moment_y", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + std::span{moment_y}); + write_numeric_1d( + context, + sections_group.get(), + "moment_z", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + std::span{moment_z}); + write_numeric_1d( + context, + sections_group.get(), + "torsion_constant", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + std::span{torsion_constants}); write_numeric_1d( context, sections_group.get(), @@ -692,7 +919,197 @@ void write_model( H5T_STD_U8LE, H5T_NATIVE_UINT8, std::span{shear_sources}); + write_numeric_2d( + context, + sections_group.get(), + "orientation", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + domain.sections().size(), + 3U, + std::span{orientations}, + "global", + "X,Y,Z"); + write_numeric_1d( + context, + sections_group.get(), + "recovery_point_offsets", + H5T_STD_U64LE, + H5T_NATIVE_UINT64, + std::span{recovery_offsets}); + write_numeric_2d( + context, + sections_group.get(), + "recovery_points", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + recovery_points.size() / 2U, + 2U, + std::span{recovery_points}, + "element_local", + "y,z"); } + + auto sets_group = create_group(context, model_group.get(), "sets"); + { + auto node_sets_group = create_group(context, sets_group.get(), "node"); + std::vector names; + std::vector offsets{0U}; + std::vector members; + names.reserve(domain.node_sets().size()); + offsets.reserve(domain.node_sets().size() + 1U); + for (const NodeSet& set : domain.node_sets()) { + names.push_back(set.name); + for (const NodeId member : set.members) { + members.push_back(member.value()); + } + offsets.push_back(static_cast(members.size())); + } + write_strings_1d(context, node_sets_group.get(), "names", names); + write_numeric_1d( + context, + node_sets_group.get(), + "member_offsets", + H5T_STD_U64LE, + H5T_NATIVE_UINT64, + std::span{offsets}); + write_numeric_1d( + context, + node_sets_group.get(), + "members", + H5T_STD_I64LE, + H5T_NATIVE_INT64, + std::span{members}); + } + { + auto element_sets_group = + create_group(context, sets_group.get(), "element"); + std::vector names; + std::vector offsets{0U}; + std::vector members; + names.reserve(domain.element_sets().size()); + offsets.reserve(domain.element_sets().size() + 1U); + for (const ElementSet& set : domain.element_sets()) { + names.push_back(set.name); + for (const ElementId member : set.members) { + members.push_back(member.value()); + } + offsets.push_back(static_cast(members.size())); + } + write_strings_1d(context, element_sets_group.get(), "names", names); + write_numeric_1d( + context, + element_sets_group.get(), + "member_offsets", + H5T_STD_U64LE, + H5T_NATIVE_UINT64, + std::span{offsets}); + write_numeric_1d( + context, + element_sets_group.get(), + "members", + H5T_STD_I64LE, + H5T_NATIVE_INT64, + std::span{members}); + } +} + +void write_analysis( + Hdf5Context& context, + const hid_t file, + const Domain& domain) { + auto analysis_group = create_group(context, file, "analysis"); + auto steps_group = create_group(context, analysis_group.get(), "steps"); + auto step_group = create_group(context, steps_group.get(), "0"); + write_string_attribute( + context, step_group.get(), "name", domain.step().name); + + auto boundary_group = + create_group(context, step_group.get(), "boundary_conditions"); + std::vector boundary_nodes; + std::vector boundary_dofs; + std::vector boundary_values; + boundary_nodes.reserve(domain.step().prescribed_dofs.size()); + boundary_dofs.reserve(domain.step().prescribed_dofs.size()); + boundary_values.reserve(domain.step().prescribed_dofs.size()); + for (const PrescribedDof& prescribed : domain.step().prescribed_dofs) { + boundary_nodes.push_back(prescribed.node.value()); + boundary_dofs.push_back(prescribed.dof); + boundary_values.push_back(prescribed.value); + } + write_numeric_1d( + context, + boundary_group.get(), + "node_ids", + H5T_STD_I64LE, + H5T_NATIVE_INT64, + std::span{boundary_nodes}); + write_numeric_1d( + context, + boundary_group.get(), + "dofs", + H5T_STD_U8LE, + H5T_NATIVE_UINT8, + std::span{boundary_dofs}); + write_numeric_1d( + context, + boundary_group.get(), + "values", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + std::span{boundary_values}); + + auto loads_group = + create_group(context, step_group.get(), "nodal_loads"); + std::vector load_nodes; + std::vector load_values; + load_nodes.reserve(domain.step().nodal_loads.size()); + load_values.reserve(domain.step().nodal_loads.size() * 6U); + for (const NodalLoad& load : domain.step().nodal_loads) { + load_nodes.push_back(load.node.value()); + load_values.insert( + load_values.end(), load.values.begin(), load.values.end()); + } + write_numeric_1d( + context, + loads_group.get(), + "node_ids", + H5T_STD_I64LE, + H5T_NATIVE_INT64, + std::span{load_nodes}); + write_numeric_2d( + context, + loads_group.get(), + "values", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + domain.step().nodal_loads.size(), + 6U, + std::span{load_values}, + "global", + "Fx,Fy,Fz,Mx,My,Mz"); + + auto solver_group = + create_group(context, analysis_group.get(), "solver_settings"); + write_string_attribute( + context, solver_group.get(), "backend", solver_backend); + write_string_attribute( + context, + solver_group.get(), + "matrix_storage", + solver_matrix_storage); + write_string_attribute( + context, + solver_group.get(), + "matrix_type", + solver_matrix_type); + write_string_attribute( + context, + solver_group.get(), + "constraint_method", + solver_constraint_method); + write_string_attribute( + context, solver_group.get(), "assembly", solver_assembly); } void write_results( @@ -752,7 +1169,9 @@ void write_results( H5T_NATIVE_DOUBLE, frame.nodal.node_ids.size(), 6U, - std::span{displacement}); + std::span{displacement}, + "global", + "Ux,Uy,Uz,Rx,Ry,Rz"); write_numeric_2d( context, nodal_group.get(), @@ -761,11 +1180,283 @@ void write_results( H5T_NATIVE_DOUBLE, frame.nodal.node_ids.size(), 6U, - std::span{reaction}); + std::span{reaction}, + "global", + "RFx,RFy,RFz,RMx,RMy,RMz"); + + auto element_group = + create_group(context, frame_group.get(), "element"); + auto beam_group = + create_group(context, element_group.get(), "beam"); + std::vector element_ids; + std::vector part_names; + std::vector instance_names; + std::vector local_labels; + std::vector local_frames; + std::vector end_node_ids; + std::vector xis; + std::vector section_strains; + std::vector section_forces; + std::vector centroid_stresses; + std::vector recovery_offsets{0U}; + std::vector recovery_stresses; + element_ids.reserve(frame.element.beams.size()); + part_names.reserve(frame.element.beams.size()); + instance_names.reserve(frame.element.beams.size()); + local_labels.reserve(frame.element.beams.size()); + local_frames.reserve(frame.element.beams.size() * 9U); + end_node_ids.reserve(frame.element.beams.size() * 2U); + xis.reserve(frame.element.beams.size() * 2U); + section_strains.reserve(frame.element.beams.size() * 12U); + section_forces.reserve(frame.element.beams.size() * 12U); + centroid_stresses.reserve(frame.element.beams.size() * 2U); + recovery_offsets.reserve(frame.element.beams.size() + 1U); + for (const BeamElementFrame& beam : frame.element.beams) { + element_ids.push_back(beam.element.value()); + part_names.push_back(beam.origin.part_name); + instance_names.push_back(beam.origin.instance_name); + local_labels.push_back(beam.origin.local_label); + for (const Vec3 axis : { + beam.local_frame.ex, + beam.local_frame.ey, + beam.local_frame.ez}) { + local_frames.push_back(axis.x); + local_frames.push_back(axis.y); + local_frames.push_back(axis.z); + } + for (const BeamSectionResult& end : beam.end_results) { + end_node_ids.push_back(end.end_node.value()); + xis.push_back(end.xi); + section_strains.insert( + section_strains.end(), + end.section_strain.begin(), + end.section_strain.end()); + section_forces.insert( + section_forces.end(), + end.section_force.begin(), + end.section_force.end()); + centroid_stresses.push_back(end.centroid_sigma_xx); + } + for (std::size_t point = 0; + point < beam.end_results[0].sigma_xx.size(); + ++point) { + recovery_stresses.push_back( + beam.end_results[0].sigma_xx[point]); + recovery_stresses.push_back( + beam.end_results[1].sigma_xx[point]); + } + recovery_offsets.push_back(static_cast( + recovery_stresses.size() / 2U)); + } + write_numeric_1d( + context, + beam_group.get(), + "element_ids", + H5T_STD_I64LE, + H5T_NATIVE_INT64, + std::span{element_ids}); + write_strings_1d( + context, beam_group.get(), "part_name", part_names); + write_strings_1d( + context, + beam_group.get(), + "instance_name", + instance_names); + write_numeric_1d( + context, + beam_group.get(), + "local_label", + H5T_STD_I64LE, + H5T_NATIVE_INT64, + std::span{local_labels}); + write_numeric_3d( + context, + beam_group.get(), + "local_frame", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + frame.element.beams.size(), + 3U, + 3U, + std::span{local_frames}, + "global", + "ex,ey,ez"); + write_numeric_2d( + context, + beam_group.get(), + "end_node_ids", + H5T_STD_I64LE, + H5T_NATIVE_INT64, + frame.element.beams.size(), + 2U, + std::span{end_node_ids}); + write_numeric_2d( + context, + beam_group.get(), + "xi", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + frame.element.beams.size(), + 2U, + std::span{xis}, + {}, + "end_minus,end_plus"); + write_numeric_3d( + context, + beam_group.get(), + "section_strain", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + frame.element.beams.size(), + 2U, + 6U, + std::span{section_strains}, + "element_local", + "epsilon,gamma_y,gamma_z,kappa_x,kappa_y,kappa_z"); + write_numeric_3d( + context, + beam_group.get(), + "section_force", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + frame.element.beams.size(), + 2U, + 6U, + std::span{section_forces}, + "element_local", + "N,Vy,Vz,T,My,Mz"); + write_numeric_2d( + context, + beam_group.get(), + "centroid_sigma_xx", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + frame.element.beams.size(), + 2U, + std::span{centroid_stresses}, + "element_local", + "end_minus,end_plus", + "sigma_xx"); + write_numeric_1d( + context, + beam_group.get(), + "recovery_point_offsets", + H5T_STD_U64LE, + H5T_NATIVE_UINT64, + std::span{recovery_offsets}); + write_numeric_2d( + context, + beam_group.get(), + "recovery_point_sigma_xx", + H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, + recovery_stresses.size() / 2U, + 2U, + std::span{recovery_stresses}, + "element_local", + "end_minus,end_plus", + "sigma_xx"); + + auto diagnostics_group = + create_group(context, frame_group.get(), "diagnostics"); + std::vector stages; + std::vector severities; + std::vector codes; + std::vector messages; + std::vector has_source; + std::vector source_files; + std::vector source_lines; + std::vector source_columns; + stages.reserve(frame.diagnostics.size()); + severities.reserve(frame.diagnostics.size()); + codes.reserve(frame.diagnostics.size()); + messages.reserve(frame.diagnostics.size()); + has_source.reserve(frame.diagnostics.size()); + source_files.reserve(frame.diagnostics.size()); + source_lines.reserve(frame.diagnostics.size()); + source_columns.reserve(frame.diagnostics.size()); + for (const Diagnostic& diagnostic : frame.diagnostics) { + stages.push_back( + static_cast(diagnostic.stage)); + severities.push_back( + static_cast(diagnostic.severity)); + codes.push_back(diagnostic.code); + messages.push_back(diagnostic.message); + has_source.push_back(diagnostic.source.has_value() ? 1U : 0U); + source_files.push_back( + diagnostic.source.has_value() + ? path_utf8(diagnostic.source->file) + : std::string{}); + source_lines.push_back( + diagnostic.source.has_value() + ? static_cast(diagnostic.source->line) + : 0U); + source_columns.push_back( + diagnostic.source.has_value() + ? static_cast( + diagnostic.source->column) + : 0U); + } + write_numeric_1d( + context, + diagnostics_group.get(), + "stage", + H5T_STD_U8LE, + H5T_NATIVE_UINT8, + std::span{stages}); + write_numeric_1d( + context, + diagnostics_group.get(), + "severity", + H5T_STD_U8LE, + H5T_NATIVE_UINT8, + std::span{severities}); + write_strings_1d( + context, diagnostics_group.get(), "code", codes); + write_strings_1d( + context, diagnostics_group.get(), "message", messages); + write_numeric_1d( + context, + diagnostics_group.get(), + "has_source", + H5T_STD_U8LE, + H5T_NATIVE_UINT8, + std::span{has_source}); + write_strings_1d( + context, + diagnostics_group.get(), + "source_file", + source_files); + write_numeric_1d( + context, + diagnostics_group.get(), + "source_line", + H5T_STD_U64LE, + H5T_NATIVE_UINT64, + std::span{source_lines}); + write_numeric_1d( + context, + diagnostics_group.get(), + "source_column", + H5T_STD_U64LE, + H5T_NATIVE_UINT64, + std::span{source_columns}); } } } +std::string read_string_attribute( + Hdf5Context& context, + hid_t parent, + std::string_view name); + +void require_field_metadata( + Hdf5Context& context, + hid_t dataset, + std::string_view coordinate_system, + std::string_view components, + std::string_view quantity = {}); + void require_numeric_type( Hdf5Context& context, const hid_t dataset, @@ -862,7 +1553,10 @@ std::vector read_numeric_2d( const std::string_view name, const hid_t file_type, const hid_t memory_type, - const std::size_t expected_columns) { + const std::size_t expected_columns, + const std::string_view coordinate_system = {}, + const std::string_view components = {}, + const std::string_view quantity = {}) { const std::string owned_name{name}; auto dataset = checked_handle( context, @@ -921,6 +1615,91 @@ std::vector read_numeric_2d( "read dataset '" + owned_name + "'", "hdf5.read_failed"); } + require_field_metadata( + context, + dataset.get(), + coordinate_system, + components, + quantity); + return values; +} + +template +std::vector read_numeric_3d( + Hdf5Context& context, + const hid_t parent, + const std::string_view name, + const hid_t file_type, + const hid_t memory_type, + const std::size_t expected_second, + const std::size_t expected_third, + const std::string_view coordinate_system = {}, + const std::string_view components = {}) { + const std::string owned_name{name}; + auto dataset = checked_handle( + context, + call_hdf5([&] { + return H5Dopen2(parent, owned_name.c_str(), H5P_DEFAULT); + }), + &H5Dclose, + "dataset '" + owned_name + "'", + "hdf5.read_failed"); + require_numeric_type(context, dataset.get(), file_type, owned_name); + auto space = checked_handle( + context, + call_hdf5([&] { return H5Dget_space(dataset.get()); }), + &H5Sclose, + "dataspace for dataset '" + owned_name + "'", + "hdf5.read_failed"); + const int rank = call_hdf5( + [&] { return H5Sget_simple_extent_ndims(space.get()); }); + if (rank != 3) { + fail( + "hdf5.read_failed", + "Dataset '" + owned_name + "' must have rank 3."); + } + std::array dimensions{}; + check_status( + call_hdf5([&] { + return H5Sget_simple_extent_dims( + space.get(), dimensions.data(), nullptr); + }), + "read shape for dataset '" + owned_name + "'", + "hdf5.read_failed"); + if ( + dimensions[1] != static_cast(expected_second) || + dimensions[2] != static_cast(expected_third)) { + fail( + "hdf5.read_failed", + "Dataset '" + owned_name + + "' has an unexpected trailing shape."); + } + const std::size_t first_count = + checked_size(dimensions[0], "Dataset row count"); + if ( + expected_second != 0U && expected_third != 0U && + first_count > std::numeric_limits::max() / + expected_second / expected_third) { + fail("hdf5.read_failed", "Dataset element count overflows."); + } + std::vector values( + first_count * expected_second * expected_third); + if (!values.empty()) { + check_status( + call_hdf5([&] { + return H5Dread( + dataset.get(), + memory_type, + H5S_ALL, + H5S_ALL, + H5P_DEFAULT, + values.data()); + }), + "read dataset '" + owned_name + "'", + "hdf5.read_failed"); + } + require_field_metadata( + context, dataset.get(), coordinate_system, components); return values; } @@ -1067,6 +1846,32 @@ std::string read_string_attribute( return std::string{raw_value}; } +void require_field_metadata( + Hdf5Context& context, + const hid_t dataset, + const std::string_view coordinate_system, + const std::string_view components, + const std::string_view quantity) { + const auto require = [&]( + const std::string_view name, + const std::string_view expected) { + if (expected.empty()) { + return; + } + const std::string actual = + read_string_attribute(context, dataset, name); + if (actual != expected) { + fail( + "hdf5.read_failed", + "Field attribute '" + std::string{name} + + "' has an unexpected value."); + } + }; + require("coordinate_system", coordinate_system); + require("components", components); + require("quantity", quantity); +} + double read_double_attribute( Hdf5Context& context, const hid_t parent, @@ -1128,6 +1933,21 @@ std::size_t group_link_count( return checked_size(information.nlinks, "Group link count"); } +void require_offsets( + const std::span offsets, + const std::size_t row_count, + const std::size_t flattened_count, + const std::string_view description) { + if ( + offsets.size() != row_count + 1U || offsets.front() != 0U || + offsets.back() != static_cast(flattened_count) || + !std::ranges::is_sorted(offsets)) { + fail( + "hdf5.read_failed", + std::string{description} + " offsets are invalid."); + } +} + Hdf5ModelSnapshot read_model( Hdf5Context& context, const hid_t file) { @@ -1135,49 +1955,33 @@ Hdf5ModelSnapshot read_model( Hdf5ModelSnapshot model; { - auto nodes_group = open_group(context, model_group.get(), "nodes"); - const auto dense_indices = read_numeric_1d( - context, - nodes_group.get(), - "dense_index", - H5T_STD_U64LE, + auto group = open_group(context, model_group.get(), "nodes"); + const auto dense = read_numeric_1d( + context, group.get(), "dense_index", H5T_STD_U64LE, H5T_NATIVE_UINT64); const auto ids = read_numeric_1d( - context, - nodes_group.get(), - "internal_id", - H5T_STD_I64LE, + context, group.get(), "internal_id", H5T_STD_I64LE, H5T_NATIVE_INT64); - const auto part_names = - read_strings_1d(context, nodes_group.get(), "part_name"); - const auto instance_names = - read_strings_1d(context, nodes_group.get(), "instance_name"); - const auto local_labels = read_numeric_1d( - context, - nodes_group.get(), - "local_label", - H5T_STD_I64LE, + const auto parts = read_strings_1d(context, group.get(), "part_name"); + const auto instances = + read_strings_1d(context, group.get(), "instance_name"); + const auto labels = read_numeric_1d( + context, group.get(), "local_label", H5T_STD_I64LE, H5T_NATIVE_INT64); const auto coordinates = read_numeric_2d( - context, - nodes_group.get(), - "coordinates", - H5T_IEEE_F64LE, - H5T_NATIVE_DOUBLE, - 3U); - const std::size_t count = dense_indices.size(); + context, group.get(), "coordinates", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, 3U, "global", "X,Y,Z"); + const std::size_t count = dense.size(); if ( - ids.size() != count || part_names.size() != count || - instance_names.size() != count || local_labels.size() != count || + ids.size() != count || parts.size() != count || + instances.size() != count || labels.size() != count || coordinates.size() != count * 3U) { - fail( - "hdf5.read_failed", - "Node datasets have inconsistent row counts."); + fail("hdf5.read_failed", "Node datasets have inconsistent rows."); } - model.nodes.reserve(count); std::unordered_set unique_ids; + model.nodes.reserve(count); for (std::size_t index = 0; index < count; ++index) { - if (dense_indices[index] != static_cast(index)) { + if (dense[index] != static_cast(index)) { fail( "hdf5.read_failed", "Node dense_index must be contiguous and 0-based."); @@ -1187,144 +1991,212 @@ Hdf5ModelSnapshot read_model( "hdf5.invalid_model_data", "Serialized model contains duplicate node internal IDs."); } - const Vec3 position{ - coordinates[index * 3U], - coordinates[index * 3U + 1U], - coordinates[index * 3U + 2U], - }; - if (!is_finite(position)) { + const Vec3 value{ + coordinates[index * 3U], coordinates[index * 3U + 1U], + coordinates[index * 3U + 2U]}; + if (!is_finite(value)) { fail( "hdf5.invalid_model_data", "Serialized model contains nonfinite node coordinates."); } model.nodes.push_back({ - dense_indices[index], - NodeId{ids[index]}, - EntityOrigin{ - part_names[index], - instance_names[index], - local_labels[index], - }, - position, - }); + dense[index], NodeId{ids[index]}, + {parts[index], instances[index], labels[index]}, value}); } } { - auto sections_group = - open_group(context, model_group.get(), "sections"); + auto group = open_group(context, model_group.get(), "materials"); const auto ids = read_numeric_1d( - context, - sections_group.get(), - "internal_id", - H5T_STD_I64LE, + context, group.get(), "internal_id", H5T_STD_I64LE, H5T_NATIVE_INT64); - const auto shear_area_y = read_numeric_1d( - context, - sections_group.get(), - "shear_area_y", - H5T_IEEE_F64LE, + const auto names = read_strings_1d(context, group.get(), "name"); + const auto young = read_numeric_1d( + context, group.get(), "young_modulus", H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE); - const auto shear_area_z = read_numeric_1d( - context, - sections_group.get(), - "shear_area_z", - H5T_IEEE_F64LE, + const auto poisson = read_numeric_1d( + context, group.get(), "poisson_ratio", H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE); - const auto shear_sources = read_numeric_1d( - context, - sections_group.get(), - "shear_source", - H5T_STD_U8LE, - H5T_NATIVE_UINT8); if ( - shear_area_y.size() != ids.size() || - shear_area_z.size() != ids.size() || - shear_sources.size() != ids.size()) { + names.size() != ids.size() || young.size() != ids.size() || + poisson.size() != ids.size()) { fail( "hdf5.read_failed", - "Section datasets have inconsistent row counts."); + "Material datasets have inconsistent rows."); } - model.sections.reserve(ids.size()); std::unordered_set unique_ids; + model.materials.reserve(ids.size()); for (std::size_t index = 0; index < ids.size(); ++index) { if (!unique_ids.insert(ids[index]).second) { fail( "hdf5.invalid_model_data", - "Serialized model contains duplicate section internal " - "IDs."); + "Serialized model contains duplicate material IDs."); } if ( - !std::isfinite(shear_area_y[index]) || - !std::isfinite(shear_area_z[index]) || - shear_area_y[index] <= 0.0 || shear_area_z[index] <= 0.0) { + !std::isfinite(young[index]) || young[index] <= 0.0 || + !std::isfinite(poisson[index]) || poisson[index] <= -1.0 || + poisson[index] >= 0.5) { fail( "hdf5.invalid_model_data", - "Serialized model requires finite positive shear areas."); + "Serialized model contains invalid material data."); + } + model.materials.push_back( + {MaterialId{ids[index]}, names[index], young[index], + poisson[index]}); + } + } + + { + auto group = open_group(context, model_group.get(), "sections"); + const auto ids = read_numeric_1d( + context, group.get(), "internal_id", H5T_STD_I64LE, + H5T_NATIVE_INT64); + const auto names = read_strings_1d(context, group.get(), "name"); + const auto areas = read_numeric_1d( + context, group.get(), "area", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE); + const auto iy = read_numeric_1d( + context, group.get(), "moment_y", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE); + const auto iz = read_numeric_1d( + context, group.get(), "moment_z", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE); + const auto torsion = read_numeric_1d( + context, group.get(), "torsion_constant", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE); + const auto shear_y = read_numeric_1d( + context, group.get(), "shear_area_y", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE); + const auto shear_z = read_numeric_1d( + context, group.get(), "shear_area_z", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE); + const auto sources = read_numeric_1d( + context, group.get(), "shear_source", H5T_STD_U8LE, + H5T_NATIVE_UINT8); + const auto orientations = read_numeric_2d( + context, group.get(), "orientation", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, 3U, "global", "X,Y,Z"); + const auto offsets = read_numeric_1d( + context, group.get(), "recovery_point_offsets", H5T_STD_U64LE, + H5T_NATIVE_UINT64); + const auto points = read_numeric_2d( + context, group.get(), "recovery_points", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, 2U, "element_local", "y,z"); + const std::size_t count = ids.size(); + if ( + names.size() != count || areas.size() != count || + iy.size() != count || iz.size() != count || + torsion.size() != count || shear_y.size() != count || + shear_z.size() != count || sources.size() != count || + orientations.size() != count * 3U) { + fail( + "hdf5.read_failed", + "Section datasets have inconsistent rows."); + } + require_offsets(offsets, count, points.size() / 2U, "Section"); + std::unordered_set unique_ids; + model.sections.reserve(count); + for (std::size_t index = 0; index < count; ++index) { + if (!unique_ids.insert(ids[index]).second) { + fail( + "hdf5.invalid_model_data", + "Serialized model contains duplicate section internal IDs."); + } + const std::array positive{ + areas[index], iy[index], iz[index], torsion[index], + shear_y[index], shear_z[index]}; + if (!std::ranges::all_of(positive, [](const double value) { + return std::isfinite(value) && value > 0.0; + })) { + fail( + "hdf5.invalid_model_data", + "Serialized model requires finite positive section data."); } ShearPropertySource source; - if (shear_sources[index] == 0U) { + if (sources[index] == 0U) { source = ShearPropertySource::input; - } else if (shear_sources[index] == 1U) { + } else if (sources[index] == 1U) { source = ShearPropertySource::phase1_default; } else { fail( "hdf5.read_failed", "Section shear_source contains an unknown value."); } + const Vec3 orientation{ + orientations[index * 3U], orientations[index * 3U + 1U], + orientations[index * 3U + 2U]}; + if (!is_finite(orientation)) { + fail( + "hdf5.invalid_model_data", + "Serialized model contains a nonfinite orientation."); + } + std::vector> recovery; + for (std::uint64_t point = offsets[index]; + point < offsets[index + 1U]; ++point) { + const std::size_t value = static_cast(point) * 2U; + if ( + !std::isfinite(points[value]) || + !std::isfinite(points[value + 1U])) { + fail( + "hdf5.invalid_model_data", + "Serialized model contains a nonfinite recovery point."); + } + recovery.push_back({points[value], points[value + 1U]}); + } model.sections.push_back({ - SectionId{ids[index]}, - shear_area_y[index], - shear_area_z[index], - source, - }); + SectionId{ids[index]}, names[index], areas[index], iy[index], + iz[index], torsion[index], shear_y[index], shear_z[index], + source, orientation, std::move(recovery)}); } } { - auto elements_group = - open_group(context, model_group.get(), "elements"); - const auto dense_indices = read_numeric_1d( - context, - elements_group.get(), - "dense_index", - H5T_STD_U64LE, + auto group = open_group(context, model_group.get(), "elements"); + const auto dense = read_numeric_1d( + context, group.get(), "dense_index", H5T_STD_U64LE, H5T_NATIVE_UINT64); const auto ids = read_numeric_1d( - context, - elements_group.get(), - "internal_id", - H5T_STD_I64LE, + context, group.get(), "internal_id", H5T_STD_I64LE, + H5T_NATIVE_INT64); + const auto parts = read_strings_1d(context, group.get(), "part_name"); + const auto instances = + read_strings_1d(context, group.get(), "instance_name"); + const auto labels = read_numeric_1d( + context, group.get(), "local_label", H5T_STD_I64LE, H5T_NATIVE_INT64); const auto connectivity = read_numeric_2d( - context, - elements_group.get(), - "connectivity", - H5T_STD_U64LE, - H5T_NATIVE_UINT64, - 2U); - const auto section_ids = read_numeric_1d( - context, - elements_group.get(), - "section_id", - H5T_STD_I64LE, + context, group.get(), "connectivity", H5T_STD_U64LE, + H5T_NATIVE_UINT64, 2U); + const auto material_ids = read_numeric_1d( + context, group.get(), "material_id", H5T_STD_I64LE, H5T_NATIVE_INT64); - const std::size_t count = dense_indices.size(); + const auto section_ids = read_numeric_1d( + context, group.get(), "section_id", H5T_STD_I64LE, + H5T_NATIVE_INT64); + const std::size_t count = dense.size(); if ( - ids.size() != count || connectivity.size() != count * 2U || - section_ids.size() != count) { + ids.size() != count || parts.size() != count || + instances.size() != count || labels.size() != count || + connectivity.size() != count * 2U || + material_ids.size() != count || section_ids.size() != count) { fail( "hdf5.read_failed", - "Element datasets have inconsistent row counts."); + "Element datasets have inconsistent rows."); } - std::unordered_set available_sections; + std::unordered_set materials; + for (const IsotropicElastic& material : model.materials) { + materials.insert(material.id.value()); + } + std::unordered_set sections; for (const Hdf5SectionSnapshot& section : model.sections) { - available_sections.insert(section.id.value()); + sections.insert(section.id.value()); } - model.elements.reserve(count); std::unordered_set unique_ids; + model.elements.reserve(count); for (std::size_t index = 0; index < count; ++index) { - if (dense_indices[index] != static_cast(index)) { + const std::uint64_t first = connectivity[index * 2U]; + const std::uint64_t second = connectivity[index * 2U + 1U]; + if (dense[index] != static_cast(index)) { fail( "hdf5.read_failed", "Element dense_index must be contiguous and 0-based."); @@ -1332,35 +2204,200 @@ Hdf5ModelSnapshot read_model( if (!unique_ids.insert(ids[index]).second) { fail( "hdf5.invalid_model_data", - "Serialized model contains duplicate element internal " - "IDs."); + "Serialized model contains duplicate element internal IDs."); } - const std::uint64_t first = connectivity[index * 2U]; - const std::uint64_t second = connectivity[index * 2U + 1U]; if ( first >= static_cast(model.nodes.size()) || - second >= static_cast(model.nodes.size())) { + second >= static_cast(model.nodes.size()) || + first == second || !materials.contains(material_ids[index]) || + !sections.contains(section_ids[index])) { fail( - "hdf5.read_failed", - "Element connectivity references a missing dense node."); - } - if (!available_sections.contains(section_ids[index])) { - fail( - "hdf5.read_failed", - "Element section_id references a missing section."); + "hdf5.invalid_model_data", + "Serialized element contains an invalid reference."); } model.elements.push_back({ - dense_indices[index], - ElementId{ids[index]}, - {first, second}, - SectionId{section_ids[index]}, - }); + dense[index], ElementId{ids[index]}, + {parts[index], instances[index], labels[index]}, + {first, second}, MaterialId{material_ids[index]}, + SectionId{section_ids[index]}}); + } + } + + auto sets_group = open_group(context, model_group.get(), "sets"); + { + auto group = open_group(context, sets_group.get(), "node"); + const auto names = read_strings_1d(context, group.get(), "names"); + const auto offsets = read_numeric_1d( + context, group.get(), "member_offsets", H5T_STD_U64LE, + H5T_NATIVE_UINT64); + const auto members = read_numeric_1d( + context, group.get(), "members", H5T_STD_I64LE, + H5T_NATIVE_INT64); + require_offsets(offsets, names.size(), members.size(), "Node set"); + std::unordered_set nodes; + for (const Hdf5NodeSnapshot& node : model.nodes) { + nodes.insert(node.id.value()); + } + for (std::size_t index = 0; index < names.size(); ++index) { + NodeSet set{names[index], {}}; + for (std::uint64_t member = offsets[index]; + member < offsets[index + 1U]; ++member) { + const std::int64_t id = + members[static_cast(member)]; + if (!nodes.contains(id)) { + fail( + "hdf5.invalid_model_data", + "Node set references a missing node."); + } + set.members.emplace_back(id); + } + model.node_sets.push_back(std::move(set)); + } + } + { + auto group = open_group(context, sets_group.get(), "element"); + const auto names = read_strings_1d(context, group.get(), "names"); + const auto offsets = read_numeric_1d( + context, group.get(), "member_offsets", H5T_STD_U64LE, + H5T_NATIVE_UINT64); + const auto members = read_numeric_1d( + context, group.get(), "members", H5T_STD_I64LE, + H5T_NATIVE_INT64); + require_offsets(offsets, names.size(), members.size(), "Element set"); + std::unordered_set elements; + for (const Hdf5ElementSnapshot& element : model.elements) { + elements.insert(element.id.value()); + } + for (std::size_t index = 0; index < names.size(); ++index) { + ElementSet set{names[index], {}}; + for (std::uint64_t member = offsets[index]; + member < offsets[index + 1U]; ++member) { + const std::int64_t id = + members[static_cast(member)]; + if (!elements.contains(id)) { + fail( + "hdf5.invalid_model_data", + "Element set references a missing element."); + } + set.members.emplace_back(id); + } + model.element_sets.push_back(std::move(set)); } } return model; } +Hdf5AnalysisSnapshot read_analysis( + Hdf5Context& context, + const hid_t file, + const Hdf5ModelSnapshot& model) { + auto analysis_group = open_group(context, file, "analysis"); + auto steps_group = open_group(context, analysis_group.get(), "steps"); + if (group_link_count(steps_group.get(), "analysis steps") != 1U) { + fail( + "hdf5.read_failed", + "Phase 1 HDF5 requires exactly one analysis step."); + } + auto step_group = open_group(context, steps_group.get(), "0"); + StepDefinition step; + step.name = read_string_attribute(context, step_group.get(), "name"); + + std::unordered_set node_ids; + for (const Hdf5NodeSnapshot& node : model.nodes) { + node_ids.insert(node.id.value()); + } + { + auto group = + open_group(context, step_group.get(), "boundary_conditions"); + const auto nodes = read_numeric_1d( + context, group.get(), "node_ids", H5T_STD_I64LE, + H5T_NATIVE_INT64); + const auto dofs = read_numeric_1d( + context, group.get(), "dofs", H5T_STD_U8LE, + H5T_NATIVE_UINT8); + const auto values = read_numeric_1d( + context, group.get(), "values", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE); + if (dofs.size() != nodes.size() || values.size() != nodes.size()) { + fail( + "hdf5.read_failed", + "Boundary-condition datasets have inconsistent rows."); + } + step.prescribed_dofs.reserve(nodes.size()); + for (std::size_t index = 0; index < nodes.size(); ++index) { + if ( + !node_ids.contains(nodes[index]) || dofs[index] < 1U || + dofs[index] > 6U || !std::isfinite(values[index])) { + fail( + "hdf5.invalid_analysis_data", + "Serialized boundary condition is invalid."); + } + step.prescribed_dofs.push_back( + {NodeId{nodes[index]}, dofs[index], values[index]}); + } + } + { + auto group = open_group(context, step_group.get(), "nodal_loads"); + const auto nodes = read_numeric_1d( + context, group.get(), "node_ids", H5T_STD_I64LE, + H5T_NATIVE_INT64); + const auto values = read_numeric_2d( + context, group.get(), "values", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, 6U, "global", "Fx,Fy,Fz,Mx,My,Mz"); + if (values.size() != nodes.size() * 6U) { + fail( + "hdf5.read_failed", + "Nodal-load datasets have inconsistent rows."); + } + step.nodal_loads.reserve(nodes.size()); + for (std::size_t index = 0; index < nodes.size(); ++index) { + if (!node_ids.contains(nodes[index])) { + fail( + "hdf5.invalid_analysis_data", + "Serialized nodal load references a missing node."); + } + std::array load{}; + std::copy_n( + values.begin() + + static_cast(index * 6U), + 6U, + load.begin()); + if (!std::ranges::all_of(load, [](const double value) { + return std::isfinite(value); + })) { + fail( + "hdf5.invalid_analysis_data", + "Serialized nodal load contains a nonfinite value."); + } + step.nodal_loads.push_back({NodeId{nodes[index]}, load}); + } + } + + auto solver_group = + open_group(context, analysis_group.get(), "solver_settings"); + Hdf5SolverSettingsSnapshot solver{ + read_string_attribute(context, solver_group.get(), "backend"), + read_string_attribute( + context, solver_group.get(), "matrix_storage"), + read_string_attribute(context, solver_group.get(), "matrix_type"), + read_string_attribute( + context, solver_group.get(), "constraint_method"), + read_string_attribute(context, solver_group.get(), "assembly"), + }; + if ( + solver.backend != solver_backend || + solver.matrix_storage != solver_matrix_storage || + solver.matrix_type != solver_matrix_type || + solver.constraint_method != solver_constraint_method || + solver.assembly != solver_assembly) { + fail( + "hdf5.invalid_analysis_data", + "Serialized solver settings are not the Phase 1 settings."); + } + return {std::move(step), std::move(solver)}; +} + ResultDatabase read_database( Hdf5Context& context, const hid_t file, @@ -1370,6 +2407,14 @@ ResultDatabase read_database( for (const Hdf5NodeSnapshot& node : model.nodes) { node_origins.emplace(node.id.value(), node.origin); } + std::unordered_map elements; + for (const Hdf5ElementSnapshot& element : model.elements) { + elements.emplace(element.id.value(), &element); + } + std::unordered_map sections; + for (const Hdf5SectionSnapshot& section : model.sections) { + sections.emplace(section.id.value(), §ion); + } ResultDatabase database{version, {}}; auto results_group = open_group(context, file, "results"); @@ -1409,14 +2454,18 @@ ResultDatabase read_database( "displacement", H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE, - 6U); + 6U, + "global", + "Ux,Uy,Uz,Rx,Ry,Rz"); const auto reaction = read_numeric_2d( context, nodal_group.get(), "reaction", H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE, - 6U); + 6U, + "global", + "RFx,RFy,RFz,RMx,RMy,RMz"); if ( displacement.size() != node_ids.size() * 6U || reaction.size() != node_ids.size() * 6U) { @@ -1450,6 +2499,219 @@ ResultDatabase read_database( frame.nodal.displacement.push_back(displacement_values); frame.nodal.reaction.push_back(reaction_values); } + + auto element_group = + open_group(context, frame_group.get(), "element"); + auto beam_group = + open_group(context, element_group.get(), "beam"); + const auto element_ids = read_numeric_1d( + context, beam_group.get(), "element_ids", H5T_STD_I64LE, + H5T_NATIVE_INT64); + const auto part_names = + read_strings_1d(context, beam_group.get(), "part_name"); + const auto instance_names = + read_strings_1d(context, beam_group.get(), "instance_name"); + const auto local_labels = read_numeric_1d( + context, beam_group.get(), "local_label", H5T_STD_I64LE, + H5T_NATIVE_INT64); + const auto local_frames = read_numeric_3d( + context, beam_group.get(), "local_frame", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, 3U, 3U, "global", "ex,ey,ez"); + const auto end_nodes = read_numeric_2d( + context, beam_group.get(), "end_node_ids", H5T_STD_I64LE, + H5T_NATIVE_INT64, 2U); + const auto xis = read_numeric_2d( + context, beam_group.get(), "xi", H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE, 2U, {}, "end_minus,end_plus"); + const auto strains = read_numeric_3d( + context, beam_group.get(), "section_strain", + H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE, 2U, 6U, + "element_local", + "epsilon,gamma_y,gamma_z,kappa_x,kappa_y,kappa_z"); + const auto forces = read_numeric_3d( + context, beam_group.get(), "section_force", + H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE, 2U, 6U, + "element_local", "N,Vy,Vz,T,My,Mz"); + const auto centroid = read_numeric_2d( + context, beam_group.get(), "centroid_sigma_xx", + H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE, 2U, "element_local", + "end_minus,end_plus", "sigma_xx"); + const auto offsets = read_numeric_1d( + context, beam_group.get(), "recovery_point_offsets", + H5T_STD_U64LE, H5T_NATIVE_UINT64); + const auto recovery = read_numeric_2d( + context, beam_group.get(), "recovery_point_sigma_xx", + H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE, 2U, "element_local", + "end_minus,end_plus", "sigma_xx"); + const std::size_t beam_count = element_ids.size(); + if ( + part_names.size() != beam_count || + instance_names.size() != beam_count || + local_labels.size() != beam_count || + local_frames.size() != beam_count * 9U || + end_nodes.size() != beam_count * 2U || + xis.size() != beam_count * 2U || + strains.size() != beam_count * 12U || + forces.size() != beam_count * 12U || + centroid.size() != beam_count * 2U) { + fail( + "hdf5.read_failed", + "Beam result datasets have inconsistent rows."); + } + require_offsets( + offsets, beam_count, recovery.size() / 2U, "Beam result"); + frame.element.beams.reserve(beam_count); + for (std::size_t index = 0; index < beam_count; ++index) { + const auto found = elements.find(element_ids[index]); + if (found == elements.end()) { + fail( + "hdf5.invalid_result_data", + "Beam result references a missing element."); + } + const Hdf5ElementSnapshot& element = *found->second; + const EntityOrigin origin{ + part_names[index], instance_names[index], + local_labels[index]}; + if (origin != element.origin) { + fail( + "hdf5.invalid_result_data", + "Beam result provenance does not match the model."); + } + const NodeId first_model = + model.nodes[static_cast( + element.connectivity[0])] + .id; + const NodeId second_model = + model.nodes[static_cast( + element.connectivity[1])] + .id; + if ( + end_nodes[index * 2U] != first_model.value() || + end_nodes[index * 2U + 1U] != second_model.value()) { + fail( + "hdf5.invalid_result_data", + "Beam result connectivity does not match the model."); + } + const auto section = sections.find(element.section.value()); + const std::size_t recovery_count = static_cast( + offsets[index + 1U] - offsets[index]); + if ( + section == sections.end() || + recovery_count != + section->second->recovery_points.size()) { + fail( + "hdf5.invalid_result_data", + "Beam recovery-point count does not match its section."); + } + const std::size_t frame_base = index * 9U; + const BeamFrame local_frame{ + {local_frames[frame_base], local_frames[frame_base + 1U], + local_frames[frame_base + 2U]}, + {local_frames[frame_base + 3U], + local_frames[frame_base + 4U], + local_frames[frame_base + 5U]}, + {local_frames[frame_base + 6U], + local_frames[frame_base + 7U], + local_frames[frame_base + 8U]}}; + std::array ends{ + BeamSectionResult{ + xis[index * 2U], NodeId{end_nodes[index * 2U]}, + {}, {}, centroid[index * 2U], {}}, + BeamSectionResult{ + xis[index * 2U + 1U], + NodeId{end_nodes[index * 2U + 1U]}, {}, {}, + centroid[index * 2U + 1U], {}}}; + for (std::size_t end = 0; end < 2U; ++end) { + const std::size_t base = (index * 2U + end) * 6U; + std::copy_n( + strains.begin() + static_cast(base), + 6U, ends[end].section_strain.begin()); + std::copy_n( + forces.begin() + static_cast(base), + 6U, ends[end].section_force.begin()); + } + for (std::uint64_t point = offsets[index]; + point < offsets[index + 1U]; ++point) { + const std::size_t base = + static_cast(point) * 2U; + ends[0].sigma_xx.push_back(recovery[base]); + ends[1].sigma_xx.push_back(recovery[base + 1U]); + } + frame.element.beams.push_back({ + ElementId{element_ids[index]}, origin, local_frame, + std::move(ends)}); + } + + auto diagnostics_group = + open_group(context, frame_group.get(), "diagnostics"); + const auto stages = read_numeric_1d( + context, diagnostics_group.get(), "stage", H5T_STD_U8LE, + H5T_NATIVE_UINT8); + const auto severities = read_numeric_1d( + context, diagnostics_group.get(), "severity", H5T_STD_U8LE, + H5T_NATIVE_UINT8); + const auto codes = + read_strings_1d(context, diagnostics_group.get(), "code"); + const auto messages = + read_strings_1d(context, diagnostics_group.get(), "message"); + const auto has_source = read_numeric_1d( + context, diagnostics_group.get(), "has_source", + H5T_STD_U8LE, H5T_NATIVE_UINT8); + const auto source_files = read_strings_1d( + context, diagnostics_group.get(), "source_file"); + const auto source_lines = read_numeric_1d( + context, diagnostics_group.get(), "source_line", + H5T_STD_U64LE, H5T_NATIVE_UINT64); + const auto source_columns = read_numeric_1d( + context, diagnostics_group.get(), "source_column", + H5T_STD_U64LE, H5T_NATIVE_UINT64); + const std::size_t diagnostic_count = stages.size(); + if ( + severities.size() != diagnostic_count || + codes.size() != diagnostic_count || + messages.size() != diagnostic_count || + has_source.size() != diagnostic_count || + source_files.size() != diagnostic_count || + source_lines.size() != diagnostic_count || + source_columns.size() != diagnostic_count) { + fail( + "hdf5.read_failed", + "Diagnostic datasets have inconsistent rows."); + } + frame.diagnostics.reserve(diagnostic_count); + for (std::size_t index = 0; + index < diagnostic_count; ++index) { + if ( + stages[index] > + static_cast( + DiagnosticStage::validation) || + severities[index] > + static_cast(Severity::error) || + has_source[index] > 1U) { + fail( + "hdf5.invalid_result_data", + "Serialized diagnostic enum value is invalid."); + } + std::optional source; + if (has_source[index] == 1U) { + source = SourceLocation{ + std::filesystem::path{source_files[index]}, + checked_size(source_lines[index], "Source line"), + checked_size( + source_columns[index], "Source column")}; + } else if ( + !source_files[index].empty() || + source_lines[index] != 0U || + source_columns[index] != 0U) { + fail( + "hdf5.invalid_result_data", + "Diagnostic without source contains source data."); + } + frame.diagnostics.push_back({ + static_cast(stages[index]), + static_cast(severities[index]), codes[index], + messages[index], std::move(source)}); + } step.frames.push_back(std::move(frame)); } database.steps.push_back(std::move(step)); @@ -1470,20 +2732,21 @@ std::vector write_hdf5( if (database.schema_version != schema_version) { return {error_diagnostic( "hdf5.unsupported_schema", - "HDF5 writer supports schema version 1.0.0 only.")}; + "HDF5 writer supports schema version 2.0.0 only.")}; } - for (const ResultStep& step : database.steps) { - for (const ResultFrame& frame : step.frames) { - if (!frame.diagnostics.empty()) { - return {error_diagnostic( - "hdf5.unsupported_result_diagnostics", - "HDF5 schema 1.0.0 does not store frame diagnostics.")}; - } - } + if ( + database.steps.size() != 1U || + database.steps[0].name != domain.step().name) { + return {error_diagnostic( + "hdf5.analysis_result_mismatch", + "HDF5 schema 2.0.0 requires one result step matching the " + "Domain step.")}; } std::unordered_set model_node_ids; + std::unordered_map model_nodes; for (const Node& node : domain.nodes()) { model_node_ids.insert(node.id.value()); + model_nodes.emplace(node.id.value(), &node); } if (const auto unknown = find_unknown_result_node(database, model_node_ids)) { @@ -1492,6 +2755,48 @@ std::vector write_hdf5( "Result node " + std::to_string(*unknown) + " is not present in the serialized model.")}; } + std::unordered_map model_elements; + for (const BeamElement& element : domain.beam_elements()) { + model_elements.emplace(element.id.value(), &element); + } + for (const ResultStep& step : database.steps) { + for (const ResultFrame& frame : step.frames) { + for (std::size_t index = 0; + index < frame.nodal.node_ids.size(); ++index) { + const auto found = + model_nodes.find(frame.nodal.node_ids[index].value()); + if ( + found == model_nodes.end() || + frame.nodal.origins[index] != found->second->origin) { + return {error_diagnostic( + "hdf5.result_provenance_mismatch", + "Nodal result provenance does not match the model.")}; + } + } + for (const BeamElementFrame& beam : frame.element.beams) { + const auto found = model_elements.find(beam.element.value()); + if (found == model_elements.end()) { + return {error_diagnostic( + "hdf5.result_element_not_in_model", + "Beam result references an element not present in " + "the serialized model.")}; + } + const BeamElement& element = *found->second; + const BeamSection& section = domain.section(element.section); + if ( + beam.origin != element.origin || + beam.end_results[0].end_node != element.nodes[0] || + beam.end_results[1].end_node != element.nodes[1] || + beam.end_results[0].sigma_xx.size() != + section.recovery_points.size()) { + return {error_diagnostic( + "hdf5.result_element_mismatch", + "Beam result provenance, connectivity, or recovery " + "points do not match the model.")}; + } + } + } + } Hdf5Context context; try { @@ -1510,7 +2815,12 @@ std::vector write_hdf5( "hdf5.create_failed"); write_string_attribute( context, file.get(), "schema_version", schema_version); + write_string_attribute( + context, file.get(), "fesa_version", version()); + write_string_attribute( + context, file.get(), "unit_policy", unit_policy); write_model(context, file.get(), domain); + write_analysis(context, file.get(), domain); write_results(context, file.get(), database); check_status( call_hdf5([&] { return H5Fflush(file.get(), H5F_SCOPE_GLOBAL); }), @@ -1546,11 +2856,30 @@ Hdf5ReadResult read_hdf5_results(const std::filesystem::path& path) { if (version != schema_version) { fail( "hdf5.unsupported_schema", - "HDF5 reader supports schema version 1.0.0 only."); + "HDF5 reader supports schema version 2.0.0 only."); + } + Hdf5MetadataSnapshot metadata{ + version, + read_string_attribute(context, file.get(), "fesa_version"), + read_string_attribute(context, file.get(), "unit_policy"), + }; + if (metadata.unit_policy != unit_policy) { + fail( + "hdf5.read_failed", + "HDF5 unit policy is not supported."); } Hdf5ModelSnapshot model = read_model(context, file.get()); + Hdf5AnalysisSnapshot analysis = + read_analysis(context, file.get(), model); ResultDatabase database = read_database(context, file.get(), version, model); + if ( + database.steps.size() != 1U || + database.steps[0].name != analysis.step.name) { + fail( + "hdf5.invalid_result_data", + "Result step does not match the analysis step."); + } const Status validation = validate_result_database(database); if (!validation.succeeded) { file.reset(); @@ -1575,6 +2904,8 @@ Hdf5ReadResult read_hdf5_results(const std::filesystem::path& path) { std::move(database), {}, std::move(model), + std::move(metadata), + std::move(analysis), }; } catch (const Hdf5Error& error) { return { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 0002328..1880b6f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -641,8 +641,14 @@ add_test( --gtest_filter=ResultRoundTrip.* ) +add_test( + NAME SelfContainedHdf5 + COMMAND "$" + --gtest_filter=SelfContainedHdf5.* +) + set_property( - TEST Hdf5 ResultRoundTrip + TEST Hdf5 ResultRoundTrip SelfContainedHdf5 PROPERTY ENVIRONMENT_MODIFICATION ${FESA_DEPENDENCY_RUNTIME_MODIFICATIONS} ) diff --git a/tests/integration/io/hdf5_results_test.cpp b/tests/integration/io/hdf5_results_test.cpp index 77347be..2e231f3 100644 --- a/tests/integration/io/hdf5_results_test.cpp +++ b/tests/integration/io/hdf5_results_test.cpp @@ -146,11 +146,36 @@ void write_double_dataset( values.data())); } +void write_root_string_attribute( + const std::filesystem::path& path, + const std::string_view name, + const std::string_view value) { + const std::string encoded_path = hdf5_path(path); + const std::string owned_name{name}; + const std::string owned_value{value}; + TestHdf5Handle file{ + H5Fopen(encoded_path.c_str(), H5F_ACC_RDWR, H5P_DEFAULT), + &H5Fclose, + }; + TestHdf5Handle attribute{ + H5Aopen(file.get(), owned_name.c_str(), H5P_DEFAULT), + &H5Aclose, + }; + TestHdf5Handle type{H5Aget_type(attribute.get()), &H5Tclose}; + const char* pointer = owned_value.c_str(); + require_hdf5_status(H5Awrite(attribute.get(), type.get(), &pointer)); +} + std::filesystem::path round_trip_path() { return std::filesystem::path{FESA_TEST_BINARY_DIR} / "Testing" / "Temporary" / "fesa-round-trip.h5"; } +std::filesystem::path self_contained_path() { + return std::filesystem::path{FESA_TEST_BINARY_DIR} / "Testing" / + "Temporary" / "fesa-self-contained.h5"; +} + fesa::Domain make_domain() { fesa::DomainBuilder builder; builder.add_node({ @@ -180,7 +205,7 @@ fesa::Domain make_domain() { 0.032, fesa::ShearPropertySource::input, fesa::Vec3{0.0, 1.0, 0.0}, - {}, + {{0.25, -0.5}, {-0.75, 0.125}}, }); builder.add_section({ fesa::SectionId{12}, @@ -209,7 +234,22 @@ fesa::Domain make_domain() { fesa::MaterialId{6}, fesa::SectionId{4}, }); - builder.set_step({"Load/Case", {}, {}}); + builder.add_node_set({"Fixed", {fesa::NodeId{42}}}); + builder.add_node_set( + {"Loaded", {fesa::NodeId{7}, fesa::NodeId{42}}}); + builder.add_element_set( + {"AllBeams", {fesa::ElementId{9}, fesa::ElementId{17}}}); + builder.set_step({ + "Load/Case", + { + {fesa::NodeId{42}, 1, 0.0}, + {fesa::NodeId{7}, 6, 0.125}, + }, + {{ + fesa::NodeId{7}, + {100.0, -200.0, 300.0, -400.0, 500.0, -600.0}, + }}, + }); auto built = std::move(builder).build(); EXPECT_TRUE(built.domain.has_value()); @@ -219,7 +259,7 @@ fesa::Domain make_domain() { fesa::ResultDatabase make_database() { return { - "1.0.0", + "2.0.0", {{ "Load/Case", {{ @@ -248,6 +288,84 @@ fesa::ResultDatabase make_database() { }; } +fesa::BeamSectionResult make_end_result( + const double xi, + const fesa::NodeId node, + const double offset) { + return { + xi, + node, + { + offset + 1.0, + offset + 2.0, + offset + 3.0, + offset + 4.0, + offset + 5.0, + offset + 6.0, + }, + { + offset + 10.0, + offset + 20.0, + offset + 30.0, + offset + 40.0, + offset + 50.0, + offset + 60.0, + }, + offset + 70.0, + {offset + 80.0, offset + 90.0}, + }; +} + +fesa::ResultDatabase make_complete_database() { + fesa::ResultDatabase database = make_database(); + auto& frame = database.steps[0].frames[0]; + frame.element.beams = { + { + fesa::ElementId{9}, + {"BeamPart", "Beam-1", 2001}, + { + {1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}, + }, + { + make_end_result(-1.0, fesa::NodeId{7}, 0.0), + make_end_result(1.0, fesa::NodeId{42}, 100.0), + }, + }, + { + fesa::ElementId{17}, + {"BeamPart", "Beam-1", 2002}, + { + {-1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, -1.0}, + }, + { + make_end_result(-1.0, fesa::NodeId{42}, 200.0), + make_end_result(1.0, fesa::NodeId{7}, 300.0), + }, + }, + }; + frame.diagnostics = { + { + fesa::DiagnosticStage::solver, + fesa::Severity::warning, + "solver.residual", + "Residual diagnostic", + fesa::SourceLocation{"beam model.inp", 41, 7}, + }, + { + fesa::DiagnosticStage::results, + fesa::Severity::error, + "results.equilibrium", + "Equilibrium diagnostic", + std::nullopt, + }, + }; + return database; +} + bool has_results_error( const std::vector& diagnostics, const std::string_view code) { @@ -276,7 +394,7 @@ TEST(ResultRoundTrip, PreservesMinimalSchemaModelAndNodalResults) { ASSERT_TRUE(read.database.has_value()); ASSERT_TRUE(read.model.has_value()); - EXPECT_EQ(read.database->schema_version, "1.0.0"); + EXPECT_EQ(read.database->schema_version, "2.0.0"); ASSERT_EQ(read.database->steps.size(), 1U); const auto& step = read.database->steps[0]; EXPECT_EQ(step.name, "Load/Case"); @@ -354,6 +472,134 @@ TEST(ResultRoundTrip, PreservesMinimalSchemaModelAndNodalResults) { fesa::ShearPropertySource::phase1_default); } +TEST(SelfContainedHdf5, PublicReaderReconstructsCompletePhase1Contract) { + const auto path = self_contained_path(); + std::filesystem::create_directories(path.parent_path()); + std::filesystem::remove(path); + + ASSERT_TRUE( + fesa::write_hdf5(path, make_domain(), make_complete_database()) + .empty()); + + const fesa::Hdf5ReadResult read = fesa::read_hdf5_results(path); + ASSERT_TRUE(read.diagnostics.empty()); + ASSERT_TRUE(read.metadata.has_value()); + ASSERT_TRUE(read.model.has_value()); + ASSERT_TRUE(read.analysis.has_value()); + ASSERT_TRUE(read.database.has_value()); + + EXPECT_EQ(read.metadata->schema_version, "2.0.0"); + EXPECT_EQ(read.metadata->fesa_version, "0.1.0"); + EXPECT_EQ( + read.metadata->unit_policy, + "consistent_input_units_no_conversion"); + + ASSERT_EQ(read.model->nodes.size(), 2U); + ASSERT_EQ(read.model->elements.size(), 2U); + EXPECT_EQ( + read.model->elements[0].origin, + (fesa::EntityOrigin{"BeamPart", "Beam-1", 2001})); + EXPECT_EQ(read.model->elements[0].material, fesa::MaterialId{6}); + ASSERT_EQ(read.model->materials.size(), 1U); + EXPECT_EQ(read.model->materials[0].name, "Steel"); + EXPECT_DOUBLE_EQ(read.model->materials[0].young, 210.0e9); + EXPECT_DOUBLE_EQ(read.model->materials[0].poisson, 0.3); + ASSERT_EQ(read.model->sections.size(), 2U); + EXPECT_EQ(read.model->sections[0].name, "InputShear"); + EXPECT_DOUBLE_EQ(read.model->sections[0].area, 0.04); + EXPECT_DOUBLE_EQ(read.model->sections[0].iy, 1.2e-4); + EXPECT_DOUBLE_EQ(read.model->sections[0].iz, 1.4e-4); + EXPECT_DOUBLE_EQ(read.model->sections[0].torsion_j, 2.0e-4); + EXPECT_DOUBLE_EQ(read.model->sections[0].orientation.y, 1.0); + EXPECT_EQ( + read.model->sections[0].recovery_points, + (std::vector>{ + {0.25, -0.5}, {-0.75, 0.125}})); + ASSERT_EQ(read.model->node_sets.size(), 2U); + EXPECT_EQ(read.model->node_sets[0].name, "Fixed"); + EXPECT_EQ( + read.model->node_sets[1].members, + (std::vector{ + fesa::NodeId{7}, fesa::NodeId{42}})); + ASSERT_EQ(read.model->element_sets.size(), 1U); + EXPECT_EQ( + read.model->element_sets[0].members, + (std::vector{ + fesa::ElementId{9}, fesa::ElementId{17}})); + + EXPECT_EQ(read.analysis->step.name, "Load/Case"); + ASSERT_EQ(read.analysis->step.prescribed_dofs.size(), 2U); + EXPECT_EQ( + read.analysis->step.prescribed_dofs[1].node, + fesa::NodeId{7}); + EXPECT_EQ(read.analysis->step.prescribed_dofs[1].dof, 6U); + EXPECT_DOUBLE_EQ( + read.analysis->step.prescribed_dofs[1].value, 0.125); + ASSERT_EQ(read.analysis->step.nodal_loads.size(), 1U); + EXPECT_DOUBLE_EQ( + read.analysis->step.nodal_loads[0].values[5], -600.0); + EXPECT_EQ(read.analysis->solver.backend, "mkl_pardiso"); + EXPECT_EQ( + read.analysis->solver.constraint_method, + "essential_dof_elimination"); + EXPECT_EQ(read.analysis->solver.assembly, "deterministic_serial"); + + const auto& frame = read.database->steps[0].frames[0]; + ASSERT_EQ(frame.element.beams.size(), 2U); + EXPECT_EQ(frame.element.beams[0].element, fesa::ElementId{9}); + EXPECT_DOUBLE_EQ(frame.element.beams[0].local_frame.ex.x, 1.0); + EXPECT_EQ( + frame.element.beams[0].end_results[0].end_node, + fesa::NodeId{7}); + EXPECT_EQ( + frame.element.beams[0].end_results[1].section_force, + (std::array{ + 110.0, 120.0, 130.0, 140.0, 150.0, 160.0})); + EXPECT_EQ( + frame.element.beams[1].end_results[0].sigma_xx, + (std::vector{280.0, 290.0})); + ASSERT_EQ(frame.diagnostics.size(), 2U); + EXPECT_EQ(frame.diagnostics[0].stage, fesa::DiagnosticStage::solver); + EXPECT_EQ(frame.diagnostics[0].severity, fesa::Severity::warning); + ASSERT_TRUE(frame.diagnostics[0].source.has_value()); + EXPECT_EQ(frame.diagnostics[0].source->file, "beam model.inp"); + EXPECT_EQ(frame.diagnostics[0].source->line, 41U); + EXPECT_EQ(frame.diagnostics[0].source->column, 7U); + EXPECT_FALSE(frame.diagnostics[1].source.has_value()); +} + +TEST(Hdf5, RejectsVersion1AfterMajorSchemaChange) { + const auto path = std::filesystem::path{FESA_TEST_BINARY_DIR} / + "Testing" / "Temporary" / "fesa-schema-1.h5"; + std::filesystem::remove(path); + auto database = make_database(); + database.schema_version = "1.0.0"; + + const auto diagnostics = + fesa::write_hdf5(path, make_domain(), database); + + EXPECT_TRUE(has_results_error(diagnostics, "hdf5.unsupported_schema")); + EXPECT_FALSE(std::filesystem::exists(path)); +} + +TEST(Hdf5, RejectsUnlistedMinorSchemaVersion) { + const auto path = std::filesystem::path{FESA_TEST_BINARY_DIR} / + "Testing" / "Temporary" / "fesa-schema-2-1.h5"; + std::filesystem::remove(path); + ASSERT_TRUE( + fesa::write_hdf5(path, make_domain(), make_database()).empty()); + write_root_string_attribute(path, "schema_version", "2.1.0"); + + const auto read = fesa::read_hdf5_results(path); + + EXPECT_FALSE(read.database.has_value()); + EXPECT_FALSE(read.model.has_value()); + EXPECT_FALSE(read.metadata.has_value()); + EXPECT_FALSE(read.analysis.has_value()); + EXPECT_TRUE(has_results_error( + read.diagnostics, "hdf5.unsupported_schema")); +} + TEST(Hdf5, RejectsInvalidResultDatabaseBeforeWriting) { const auto path = std::filesystem::path{FESA_TEST_BINARY_DIR} / "Testing" / "Temporary" / "fesa-invalid.h5"; @@ -369,7 +615,7 @@ TEST(Hdf5, RejectsInvalidResultDatabaseBeforeWriting) { EXPECT_FALSE(std::filesystem::exists(path)); } -TEST(Hdf5, RejectsFrameDiagnosticsThatSchemaCannotRepresent) { +TEST(Hdf5, PreservesFrameDiagnosticsRepresentedBySchema) { const auto path = std::filesystem::path{FESA_TEST_BINARY_DIR} / "Testing" / "Temporary" / "fesa-unsupported-diagnostics.h5"; @@ -386,9 +632,18 @@ TEST(Hdf5, RejectsFrameDiagnosticsThatSchemaCannotRepresent) { const auto diagnostics = fesa::write_hdf5(path, domain, database); - EXPECT_TRUE(has_results_error( - diagnostics, "hdf5.unsupported_result_diagnostics")); - EXPECT_FALSE(std::filesystem::exists(path)); + ASSERT_TRUE(diagnostics.empty()); + + const auto read = fesa::read_hdf5_results(path); + ASSERT_TRUE(read.diagnostics.empty()); + ASSERT_TRUE(read.database.has_value()); + const auto& stored = read.database->steps[0].frames[0].diagnostics; + ASSERT_EQ(stored.size(), 1U); + EXPECT_EQ(stored[0].stage, fesa::DiagnosticStage::solver); + EXPECT_EQ(stored[0].severity, fesa::Severity::warning); + EXPECT_EQ(stored[0].code, "solver.residual"); + EXPECT_EQ(stored[0].message, "Residual diagnostic"); + EXPECT_FALSE(stored[0].source.has_value()); } TEST(Hdf5, RejectsResultNodeMissingFromDomainBeforeWriting) { diff --git a/tests/integration/pipeline/minimal_cantilever_test.cpp b/tests/integration/pipeline/minimal_cantilever_test.cpp index 913892e..3d6f6b4 100644 --- a/tests/integration/pipeline/minimal_cantilever_test.cpp +++ b/tests/integration/pipeline/minimal_cantilever_test.cpp @@ -94,7 +94,7 @@ TEST(MinimalCantileverPipeline, WritesReadableFiniteEquilibratedResults) { ASSERT_TRUE(read.database.has_value()); ASSERT_TRUE(read.model.has_value()); EXPECT_TRUE(read.diagnostics.empty()); - EXPECT_EQ(read.database->schema_version, "1.0.0"); + EXPECT_EQ(read.database->schema_version, "2.0.0"); ASSERT_EQ(read.model->nodes.size(), 2U); EXPECT_EQ(read.model->nodes[0].id, fesa::NodeId{0}); diff --git a/tests/unit/analysis/linear_static_analysis_test.cpp b/tests/unit/analysis/linear_static_analysis_test.cpp index 9746c7c..a971de7 100644 --- a/tests/unit/analysis/linear_static_analysis_test.cpp +++ b/tests/unit/analysis/linear_static_analysis_test.cpp @@ -149,7 +149,7 @@ TEST(LinearStaticAnalysis, SolvesHandCalculatedAxialBeamInNodeIdOrder) { ASSERT_TRUE(run.results.has_value()); EXPECT_TRUE(run.diagnostics.empty()); EXPECT_TRUE(fesa::validate_result_database(*run.results).succeeded); - EXPECT_EQ(run.results->schema_version, "1.0.0"); + EXPECT_EQ(run.results->schema_version, "2.0.0"); ASSERT_EQ(run.results->steps.size(), 1); EXPECT_EQ(run.results->steps[0].name, "Load"); ASSERT_EQ(run.results->steps[0].frames.size(), 1); diff --git a/tests/unit/results/result_database_test.cpp b/tests/unit/results/result_database_test.cpp index 74cb29f..30f5de3 100644 --- a/tests/unit/results/result_database_test.cpp +++ b/tests/unit/results/result_database_test.cpp @@ -64,7 +64,7 @@ fesa::ResultDatabase valid_database() { {}, }; fesa::ResultStep step{"Load", {std::move(frame)}}; - return {"1.0.0", {std::move(step)}}; + return {"2.0.0", {std::move(step)}}; } bool has_diagnostic(