diff --git a/include/fesa/analysis/newton_iteration_log.h b/include/fesa/analysis/newton_iteration_log.h index e2de1da..c59ab27 100644 --- a/include/fesa/analysis/newton_iteration_log.h +++ b/include/fesa/analysis/newton_iteration_log.h @@ -11,6 +11,8 @@ #include "fesa/core/status.h" namespace fesa { +struct SolidStaticIterationRecord; +struct SolidDynamicIterationRecord; /// @brief Identifies the latest Newton procedure stage reached by a record. enum class NewtonIterationStage { @@ -93,6 +95,12 @@ class NewtonIterationLogWriter { /// evidence or stream failure permanently fails this writer instance. Status Write(const DynamicNewtonIterationLogRecord& record); + /// @brief Writes the two actual Solid static convergence channels. + Status Write(const SolidStaticIterationRecord& record, + const std::string& step_name); + /// @brief Writes scalar full-pair Solid HHT convergence evidence. + Status Write(const SolidDynamicIterationRecord& record); + /// @brief Flushes every complete record written so far to the log. /// @return Success after the stream flush, otherwise an output failure with /// code `log-output-failure`. diff --git a/include/fesa/analysis/prepared_analysis.h b/include/fesa/analysis/prepared_analysis.h index 1c7a7bc..becbd3d 100644 --- a/include/fesa/analysis/prepared_analysis.h +++ b/include/fesa/analysis/prepared_analysis.h @@ -9,6 +9,9 @@ #include "fesa/analysis/dynamic_analysis_state.h" #include "fesa/analysis/eigen_analysis_state.h" #include "fesa/analysis/nonlinear_analysis_state.h" +#include "fesa/analysis/solid_frequency_state.h" +#include "fesa/analysis/solid_hht_dynamic_procedure.h" +#include "fesa/analysis/solid_static_procedure.h" #include "fesa/analysis/transient_heat_procedure.h" #include "fesa/core/diagnostic.h" #include "fesa/core/status.h" @@ -26,6 +29,8 @@ class ResultsWriter; class HhtDynamicProcedureExecutor; class DynamicResultsWriter; class ThermalResultsWriter; +class SolidResultsWriter; +class SolidFrequencyProcedureExecutor; struct DerivedOutputPaths; /// @brief Owns one parsed and semantically mapped analysis input. @@ -63,7 +68,8 @@ struct DynamicExecutionResult { using PreparedAnalysisOutcome = std::variant; + ThermalExecutionResult, SolidStaticExecutionResult, + SolidFrequencyState, SolidDynamicExecutionResult>; /// @brief Selects the approved procedure from a prepared Domain. /// @pre The executors and log outlive this call. @@ -76,7 +82,10 @@ Result ExecutePreparedAnalysis( const LinearFrequencyProcedureExecutor& frequency_executor, NewtonIterationLogWriter& log, const HhtDynamicProcedureExecutor* dynamic_executor = nullptr, - const TransientHeatProcedureExecutor* thermal_executor = nullptr); + const TransientHeatProcedureExecutor* thermal_executor = nullptr, + const SolidStaticProcedureExecutor* solid_static_executor = nullptr, + const SolidFrequencyProcedureExecutor* solid_frequency_executor = nullptr, + const SolidHhtDynamicProcedureExecutor* solid_dynamic_executor = nullptr); /// @brief Writes an already computed outcome to its authoritative result file. /// @pre The Newton log has closed successfully before this function is called. @@ -89,7 +98,8 @@ Status FinalizePreparedAnalysisResult( const std::vector& diagnostics, ResultsWriter& linear_writer, NonlinearResultsWriter& nonlinear_writer, EigenResultsWriter& eigen_writer, DynamicResultsWriter* dynamic_writer = nullptr, - ThermalResultsWriter* thermal_writer = nullptr); + ThermalResultsWriter* thermal_writer = nullptr, + SolidResultsWriter* solid_writer = nullptr); } // namespace fesa diff --git a/include/fesa/io/hdf5/solid_hdf5_results_writer.h b/include/fesa/io/hdf5/solid_hdf5_results_writer.h new file mode 100644 index 0000000..299cbc2 --- /dev/null +++ b/include/fesa/io/hdf5/solid_hdf5_results_writer.h @@ -0,0 +1,21 @@ +#ifndef FESA_IO_HDF5_SOLID_HDF5_RESULTS_WRITER_H_ +#define FESA_IO_HDF5_SOLID_HDF5_RESULTS_WRITER_H_ +#include "fesa/results/solid_results_writer.h" +namespace fesa { +/// @brief Writes and validates complete schema1 Solid RN64 outcomes atomically. +/// @note Invalid payloads and write, readback, close or replacement failures +/// return an output failure while preserving any existing final result. +class SolidHdf5ResultsWriter final : public SolidResultsWriter { + public: + Status Write(const std::filesystem::path&, const Domain&, + const SolidStaticState&, + const std::vector&) override; + Status Write(const std::filesystem::path&, const Domain&, + const SolidFrequencyState&, + const std::vector&) override; + Status Write(const std::filesystem::path&, const Domain&, + const SolidDynamicAnalysisState&, + const std::vector&) override; +}; +} // namespace fesa +#endif diff --git a/include/fesa/results/solid_results_writer.h b/include/fesa/results/solid_results_writer.h new file mode 100644 index 0000000..45039fd --- /dev/null +++ b/include/fesa/results/solid_results_writer.h @@ -0,0 +1,33 @@ +#ifndef FESA_RESULTS_SOLID_RESULTS_WRITER_H_ +#define FESA_RESULTS_SOLID_RESULTS_WRITER_H_ +#include +#include + +#include "fesa/analysis/solid_dynamic_analysis_state.h" +#include "fesa/analysis/solid_frequency_state.h" +#include "fesa/analysis/solid_static_state.h" +#include "fesa/core/diagnostic.h" +#include "fesa/core/status.h" +#include "fesa/model/domain.h" +namespace fesa { +/// @brief Persists only complete typed Solid outcomes, preserving previous +/// finals. +class SolidResultsWriter { + public: + virtual ~SolidResultsWriter() = default; + /// @brief Writes a complete static history in stable source order. + virtual Status Write(const std::filesystem::path&, const Domain&, + const SolidStaticState&, + const std::vector&) = 0; + /// @brief Writes validated modal vectors and all returned-pair evidence. + virtual Status Write(const std::filesystem::path&, const Domain&, + const SolidFrequencyState&, + const std::vector&) = 0; + /// @brief Writes nearest Float64 observables of complete physical pairs. + /// @note No recovered quantity is recomputed from rounded displacement. + virtual Status Write(const std::filesystem::path&, const Domain&, + const SolidDynamicAnalysisState&, + const std::vector&) = 0; +}; +} // namespace fesa +#endif diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index 58d1627..b8c889e 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -71,6 +71,8 @@ add_library( io/hdf5/hdf5_self_check.cpp io/hdf5/thermal_hdf5_results_writer.cpp io/hdf5/thermal_hdf5_self_check.cpp + io/hdf5/solid_hdf5_results_writer.cpp + io/hdf5/solid_hdf5_self_check.cpp loads/concentrated_nodal_load.cpp loads/shell_surface_pressure.cpp materials/isotropic_linear_elastic_material.cpp diff --git a/src/fesa/analysis/newton_iteration_log.cpp b/src/fesa/analysis/newton_iteration_log.cpp index 32fb67f..895364e 100644 --- a/src/fesa/analysis/newton_iteration_log.cpp +++ b/src/fesa/analysis/newton_iteration_log.cpp @@ -9,6 +9,9 @@ #include #include +#include "fesa/analysis/solid_dynamic_analysis_state.h" +#include "fesa/analysis/solid_static_state.h" + namespace fesa { namespace { @@ -320,6 +323,119 @@ Result SerializeRecord( return Result::Success(std::move(line)); } +Result JoinSolidFields(const std::vector& fields) { + std::string line; + for (const auto& field : fields) { + if (!line.empty()) line += '\t'; + line += field; + } + line += '\n'; + return Result::Success(std::move(line)); +} +Result SerializeRecord(const SolidStaticIterationRecord& record, + const std::string& step, + const std::filesystem::path& path) { + std::vector fields; + AppendField(&fields, "schema", "solid-static-log-v1"); + AppendField(&fields, "step", EscapeText(step)); + AppendField(&fields, "increment", std::to_string(record.increment_index)); + AppendField(&fields, "iteration", std::to_string(record.iteration)); + const auto factor = FormatFiniteNumber(record.load_factor, path); + if (!factor.HasValue()) + return Result::Failure(factor.GetStatus()); + AppendField(&fields, "load_factor", factor.Value()); + AppendField(&fields, "stage", EscapeText(record.stage)); + for (const auto& entry : + std::array, 2>{ + {{{"residual_force", "residual_force"}, &record.force_residual}, + {{"correction_translation", "correction_translation"}, + &record.translation_correction}}}) { + const Status status = AppendChannelFields( + &fields, entry.first, *entry.second, record.convergence_evaluated, + "not-evaluated(" + EscapeText(record.failure_reason) + ")", path); + if (!status.IsOk()) return Result::Failure(status); + } + AppendField(&fields, "accepted", record.accepted ? "1" : "0"); + AppendField(&fields, "diagnostic_code", EscapeText(record.diagnostic_code)); + AppendField(&fields, "failure_reason", EscapeText(record.failure_reason)); + return JoinSolidFields(fields); +} +Result SerializeRecord(const SolidDynamicIterationRecord& record, + const std::filesystem::path& path) { + std::vector fields; + AppendField(&fields, "schema", "solid-hht-log-v1"); + AppendField(&fields, "arithmetic", + "normalized-high-low-physical-state;RN64-observables"); + AppendField(&fields, "increment", std::to_string(record.increment_index)); + AppendField(&fields, "iteration", + std::to_string(record.evidence.iteration_index)); + AppendField(&fields, "stage", EscapeText(record.stage)); + AppendField(&fields, "outcome", EscapeText(record.outcome)); + AppendField(&fields, "failure_code", EscapeText(record.failure_code)); + const auto& e = record.evidence; + for (const double metric : + {e.residual_numerator, e.residual_scale, e.residual_metric, + e.weighted_residual_numerator, e.recombined_residual_numerator, + e.old_residual_numerator, e.weighted_equivalence_metric, + e.update_metric}) { + if (!std::isfinite(metric) || metric < 0) + return Result::Failure(LogOutputFailure( + path, "A Solid HHT norm, scale or metric is invalid.")); + } + for (const double metric : e.correction_metrics) { + if (!std::isfinite(metric) || metric < 0) + return Result::Failure( + LogOutputFailure(path, "A Solid HHT correction metric is invalid.")); + } + std::vector> numbers{ + {"time", record.time}, + {"physical_residual_force_norm", e.residual_numerator}, + {"uncancelled_force_scale", e.residual_scale}, + {"physical_residual_max_ratio", e.residual_metric}, + {"original_weighted_force_norm", e.weighted_residual_numerator}, + {"solved_weighted_force_norm", e.recombined_residual_numerator}, + {"old_physical_force_norm", e.old_residual_numerator}, + {"weighted_equivalence_metric", e.weighted_equivalence_metric}, + {"update_metric", e.update_metric}}; + const std::array corrections{{"du_over_L", "h_dv_over_L", + "h2_physical_da_over_L", + "h2_algorithmic_dA_over_L"}}; + for (std::size_t i = 0; i < 4; ++i) + numbers.emplace_back(corrections[i], e.correction_metrics[i]); + const std::array terms{ + {"worst_r_force", "worst_Ma_force", "worst_internal_force", + "worst_external_force", "worst_uncancelled_force_scale"}}; + for (std::size_t i = 0; i < 5; ++i) + numbers.emplace_back(terms[i], e.worst_force_terms[i]); + for (const auto& number : numbers) { + const auto formatted = FormatFiniteNumber(number.second, path); + if (!formatted.HasValue()) + return Result::Failure(formatted.GetStatus()); + AppendField(&fields, number.first, formatted.Value()); + } + const std::array weighted{ + {"worst_original_weighted_force", "worst_solved_weighted_force", + "worst_old_physical_force", "worst_weighted_uncancelled_force_scale"}}; + for (std::size_t i = 0; i < 4; ++i) { + const auto formatted = FormatFiniteNumber(e.worst_weighted_terms[i], path); + if (!formatted.HasValue()) + return Result::Failure(formatted.GetStatus()); + AppendField(&fields, weighted[i], formatted.Value()); + } + AppendField(&fields, "physical_residual_threshold", "1e-10"); + AppendField(&fields, "correction_threshold", "1e-10"); + AppendField(&fields, "update_threshold", "1e-12"); + AppendField(&fields, "weighted_equivalence_threshold", "1e-12"); + AppendField(&fields, "residual_unit", "force"); + AppendField(&fields, "correction_unit", "dimensionless"); + AppendField(&fields, "worst_full_dof", std::to_string(e.worst_full_dof)); + AppendField(&fields, "zero_scale_exact", e.zero_scale_exact ? "1" : "0"); + AppendField(&fields, "physical_residual_pass", e.residual_pass ? "1" : "0"); + AppendField(&fields, "correction_pass", e.correction_pass ? "1" : "0"); + AppendField(&fields, "update_pass", e.update_pass ? "1" : "0"); + return JoinSolidFields(fields); +} + } // namespace Status NewtonIterationLogWriter::Open(const std::filesystem::path& log_path) { @@ -379,6 +495,48 @@ Status NewtonIterationLogWriter::Write( return Status::Ok(); } +Status NewtonIterationLogWriter::Write(const SolidStaticIterationRecord& record, + const std::string& step_name) { + if (state_ != State::kOpen) { + state_ = State::kFailed; + return LogOutputFailure(log_path_, "The Newton log writer is not open."); + } + const auto line = SerializeRecord(record, step_name, log_path_); + if (!line.HasValue()) { + state_ = State::kFailed; + return line.GetStatus(); + } + stream_.write(line.Value().data(), + static_cast(line.Value().size())); + if (!stream_) { + state_ = State::kFailed; + return LogOutputFailure(log_path_, + "The Solid Newton log record could not write."); + } + return Status::Ok(); +} + +Status NewtonIterationLogWriter::Write( + const SolidDynamicIterationRecord& record) { + if (state_ != State::kOpen) { + state_ = State::kFailed; + return LogOutputFailure(log_path_, "The Newton log writer is not open."); + } + const auto line = SerializeRecord(record, log_path_); + if (!line.HasValue()) { + state_ = State::kFailed; + return line.GetStatus(); + } + stream_.write(line.Value().data(), + static_cast(line.Value().size())); + if (!stream_) { + state_ = State::kFailed; + return LogOutputFailure(log_path_, + "The Solid Newton log record could not write."); + } + return Status::Ok(); +} + Status NewtonIterationLogWriter::Flush() { if (state_ != State::kOpen) { state_ = State::kFailed; diff --git a/src/fesa/analysis/prepared_analysis.cpp b/src/fesa/analysis/prepared_analysis.cpp index 7ea59d5..fff5a2e 100644 --- a/src/fesa/analysis/prepared_analysis.cpp +++ b/src/fesa/analysis/prepared_analysis.cpp @@ -2,11 +2,13 @@ #include +#include "fesa/analysis/analysis_model.h" #include "fesa/analysis/hht_dynamic_procedure.h" #include "fesa/analysis/linear_frequency_procedure.h" #include "fesa/analysis/linear_static_procedure.h" #include "fesa/analysis/newton_iteration_log.h" #include "fesa/analysis/nonlinear_static_procedure.h" +#include "fesa/analysis/solid_frequency_procedure.h" #include "fesa/app/output_paths.h" #include "fesa/io/abaqus/domain_mapper.h" #include "fesa/io/abaqus/input_reader.h" @@ -14,6 +16,7 @@ #include "fesa/results/eigen_results_writer.h" #include "fesa/results/nonlinear_results_writer.h" #include "fesa/results/results_writer.h" +#include "fesa/results/solid_results_writer.h" #include "fesa/results/thermal_results_writer.h" namespace fesa { @@ -42,7 +45,10 @@ Result ExecutePreparedAnalysis( const LinearFrequencyProcedureExecutor& frequency_executor, NewtonIterationLogWriter& log, const HhtDynamicProcedureExecutor* dynamic_executor, - const TransientHeatProcedureExecutor* thermal_executor) { + const TransientHeatProcedureExecutor* thermal_executor, + const SolidStaticProcedureExecutor* solid_static_executor, + const SolidFrequencyProcedureExecutor* solid_frequency_executor, + const SolidHhtDynamicProcedureExecutor* solid_dynamic_executor) { if (domain.ThermalStep() != nullptr) { if (thermal_executor == nullptr) { return Result::Failure( @@ -72,6 +78,63 @@ Result ExecutePreparedAnalysis( "Prepared analysis requires exactly one analysis step."}})); } + if (domain.IsSolid()) { + const auto missing = []() { + return Result::Failure(Status::Failure( + FailureCategory::kSolver, + {{Severity::kError, + "missing-solid-executor", + {}, + "", + {}, + "The selected Solid execution adapter is required."}})); + }; + if (domain.Dynamic() != nullptr) { + if (solid_dynamic_executor == nullptr) return missing(); + auto model = AnalysisModel::Create(domain); + if (!model.HasValue()) + return Result::Failure(model.GetStatus()); + auto dofs = DofManager::Create(model.Value()); + if (!dofs.HasValue()) + return Result::Failure(dofs.GetStatus()); + auto state = SolidDynamicAnalysisState::Create(domain, dofs.Value()); + if (!state.HasValue()) + return Result::Failure(state.GetStatus()); + const Status execution = + solid_dynamic_executor->ExecuteInto(domain, state.Value()); + for (const auto& record : state.Value().IterationRecords()) { + const Status written = log.Write(record); + if (!written.IsOk()) + return Result::Failure(written); + } + if (!execution.IsOk()) + return Result::Failure(execution); + return Result::Success( + PreparedAnalysisOutcome{SolidDynamicExecutionResult{ + std::move(dofs.Value()), std::move(state.Value())}}); + } + if (domain.Steps()[0U].IsFrequency()) { + if (solid_frequency_executor == nullptr) return missing(); + auto result = solid_frequency_executor->Execute(domain); + if (!result.HasValue()) + return Result::Failure(result.GetStatus()); + return Result::Success( + PreparedAnalysisOutcome{std::move(result.Value())}); + } + if (solid_static_executor == nullptr) return missing(); + SolidStaticState state; + const Status execution = solid_static_executor->ExecuteInto(domain, state); + for (const auto& record : state.IterationRecords()) { + const Status written = log.Write(record, domain.Steps()[0U].Name()); + if (!written.IsOk()) + return Result::Failure(written); + } + if (!execution.IsOk()) + return Result::Failure(execution); + return Result::Success(PreparedAnalysisOutcome{ + SolidStaticExecutionResult{std::move(state), domain.Warnings()}}); + } + if (domain.Dynamic() != nullptr) { if (dynamic_executor == nullptr) return Result::Failure( @@ -120,8 +183,29 @@ Status FinalizePreparedAnalysisResult( const PreparedAnalysisOutcome& outcome, const std::vector& diagnostics, ResultsWriter& linear_writer, NonlinearResultsWriter& nonlinear_writer, EigenResultsWriter& eigen_writer, - DynamicResultsWriter* dynamic_writer, - ThermalResultsWriter* thermal_writer) { + DynamicResultsWriter* dynamic_writer, ThermalResultsWriter* thermal_writer, + SolidResultsWriter* solid_writer) { + if (std::holds_alternative(outcome) || + std::holds_alternative(outcome) || + std::holds_alternative(outcome)) { + if (solid_writer == nullptr) + return Status::Failure(FailureCategory::kOutput, + {{Severity::kError, + "missing-solid-writer", + {}, + "", + {}, + "Solid result writer is required."}}); + if (const auto* result = std::get_if(&outcome)) + return solid_writer->Write(paths.results_path, domain, result->state, + diagnostics); + if (const auto* result = std::get_if(&outcome)) + return solid_writer->Write(paths.results_path, domain, *result, + diagnostics); + return solid_writer->Write( + paths.results_path, domain, + std::get(outcome).state, diagnostics); + } if (std::holds_alternative(outcome)) { if (thermal_writer == nullptr) { return Status::Failure(FailureCategory::kOutput, diff --git a/src/fesa/app/fesa_application.cpp b/src/fesa/app/fesa_application.cpp index 5532552..f160877 100644 --- a/src/fesa/app/fesa_application.cpp +++ b/src/fesa/app/fesa_application.cpp @@ -11,6 +11,7 @@ #include "fesa/analysis/newton_iteration_log.h" #include "fesa/analysis/nonlinear_static_procedure.h" #include "fesa/analysis/prepared_analysis.h" +#include "fesa/analysis/solid_frequency_procedure.h" #include "fesa/app/output_paths.h" #include "fesa/assembly/parallel_for.h" #include "fesa/core/diagnostic.h" @@ -18,6 +19,7 @@ #include "fesa/io/hdf5/eigen_hdf5_results_writer.h" #include "fesa/io/hdf5/hdf5_results_writer.h" #include "fesa/io/hdf5/nonlinear_hdf5_results_writer.h" +#include "fesa/io/hdf5/solid_hdf5_results_writer.h" #include "fesa/io/hdf5/thermal_hdf5_results_writer.h" #include "fesa/solvers/linear/mkl_pardiso_solver.h" #include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" @@ -89,6 +91,7 @@ int FesaApplication::Run(const std::vector& arguments) { EigenHdf5ResultsWriter eigen_writer; DynamicHdf5ResultsWriter dynamic_writer; ThermalHdf5ResultsWriter thermal_writer; + SolidHdf5ResultsWriter solid_writer; NewtonIterationLogWriter log; Status status = log.Open(paths.Value().log_path); @@ -110,9 +113,15 @@ int FesaApplication::Run(const std::vector& arguments) { parallel_for, linear_solver, tangent_solver); const TransientHeatProcedureExecutor thermal_executor(parallel_for, linear_solver); + const SolidStaticProcedureExecutor solid_static(parallel_for, linear_solver, + tangent_solver); + const SolidFrequencyProcedureExecutor solid_frequency; + const SolidHhtDynamicProcedureExecutor solid_dynamic( + parallel_for, linear_solver, tangent_solver); Result outcome = ExecutePreparedAnalysis( prepared.Value().domain, linear_executor, nonlinear_executor, - frequency_executor, log, &dynamic_executor, &thermal_executor); + frequency_executor, log, &dynamic_executor, &thermal_executor, + &solid_static, &solid_frequency, &solid_dynamic); if (!outcome.HasValue()) { status = outcome.GetStatus(); } else { @@ -124,8 +133,10 @@ int FesaApplication::Run(const std::vector& arguments) { status = FinalizePreparedAnalysisResult( paths.Value(), prepared.Value().domain, outcome.Value(), prepared.Value().diagnostics, linear_writer, nonlinear_writer, - eigen_writer, &dynamic_writer, &thermal_writer); + eigen_writer, &dynamic_writer, &thermal_writer, &solid_writer); if (status.IsOk()) { + if (prepared.Value().domain.IsSolid()) + WriteDiagnostics(prepared.Value().diagnostics); return kSuccessExitCode; } WriteDiagnostics(status.Diagnostics()); diff --git a/src/fesa/io/hdf5/solid_hdf5_results_writer.cpp b/src/fesa/io/hdf5/solid_hdf5_results_writer.cpp new file mode 100644 index 0000000..bb6772b --- /dev/null +++ b/src/fesa/io/hdf5/solid_hdf5_results_writer.cpp @@ -0,0 +1,763 @@ +#include "fesa/io/hdf5/solid_hdf5_results_writer.h" + +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/build_info.h" +#include "fesa/elements/element_factory.h" +#include "fesa/elements/hex8_element.h" +#include "fesa/materials/isotropic_linear_elastic_material.h" +#include "io/hdf5/hdf5_atomic_file.h" +#include "io/hdf5/hdf5_model_writer.h" +#include "io/hdf5/hdf5_primitives.h" +#include "io/hdf5/hdf5_result_writer.h" +#include "io/hdf5/solid_hdf5_self_check.h" +#include "math/high_low_arithmetic.h" + +namespace fesa::hdf5_internal { +namespace { +void Check(bool condition, const char* message) { + if (!condition) throw Hdf5Failure(message); +} +using U = std::vector; +using B = std::vector; +using D = std::vector; +using S = std::vector; +std::vector Values(const Vector& vector) { + D result; + for (std::size_t i = 0; i < vector.Size(); ++i) result.push_back(vector[i]); + return result; +} +void Parent(hid_t file, const std::string& path) { + const auto parent = path.substr(0, path.find_last_of('/')); + if (!parent.empty() && H5Lexists(file, parent.c_str(), H5P_DEFAULT) <= 0) + (void)CreateGroup(file, parent.c_str()); +} +void Table(hid_t file, const SolidTable& table) { + Parent(file, table.path); + std::vector native, filetypes; + std::size_t size = 0; + for (const auto& column : table.columns) { + native.push_back(SolidColumnType(column, true)); + filetypes.push_back(SolidColumnType(column, false)); + size += H5Tget_size(native.back().Get()); + } + Hdf5Handle ft{RequireHdf5Id(H5Tcreate(H5T_COMPOUND, size), + "Solid file compound failed."), + H5Tclose}; + Hdf5Handle mt{RequireHdf5Id(H5Tcreate(H5T_COMPOUND, size), + "Solid memory compound failed."), + H5Tclose}; + std::vector bytes(table.rows * size); + std::size_t offset = 0; + for (std::size_t c = 0; c < table.columns.size(); ++c) { + const auto& column = table.columns[c]; + RequireHdf5( + H5Tinsert(ft.Get(), column.name.c_str(), offset, filetypes[c].Get()), + "Solid file member failed."); + RequireHdf5( + H5Tinsert(mt.Get(), column.name.c_str(), offset, native[c].Get()), + "Solid memory member failed."); + std::visit( + [&](const auto& values) { + using T = typename std::decay_t::value_type; + for (std::size_t row = 0; row < table.rows; ++row) { + if constexpr (std::is_same_v) { + const char* pointer = values[row].c_str(); + std::memcpy(bytes.data() + row * size + offset, &pointer, + sizeof(pointer)); + } else + std::memcpy(bytes.data() + row * size + offset, + values.data() + row * column.width, + column.width * sizeof(T)); + } + }, + column.values); + offset += H5Tget_size(native[c].Get()); + } + auto dataset = WriteCompoundDataset(file, table.path.c_str(), table.rows, + ft.Get(), mt.Get(), bytes.data()); + for (const auto& item : table.strings) + WriteStringAttribute(dataset.Get(), item.first.c_str(), item.second); +} +SolidData Data(const std::string& path, std::vector dimensions, + D values, const char* names, const char* units, + const char* location) { + return {path, + std::move(dimensions), + std::move(values), + {{"component_names", names}, + {"unit_dimensions", units}, + {"coordinate_system", "global-cartesian"}, + {"location", location}, + {"step_name", "Step-1"}}, + {}, + {}}; +} +void FrameData(SolidPlan& plan, const Domain& domain, std::size_t frame, + std::size_t increment, double coordinate, bool dynamic, + const Vector& u, const Vector& rf, const Vector& r, + const std::vector& points, + double energy) { + const std::string root = "/steps/Step-1/frames/" + std::to_string(frame); + const hsize_t n = domain.Nodes().size(), e = domain.Elements().Size(); + const auto add = [&](SolidData value) { + value.integers = {{"frame_index", frame}, {"increment_index", increment}}; + value.doubles = {{dynamic ? "time" : "load_factor", coordinate}}; + plan.datasets.push_back(std::move(value)); + }; + add(Data(root + "/nodal/displacement", {n, 3}, Values(u), "U1,U2,U3", + "length", "nodal")); + add(Data(root + "/nodal/reaction", {n, 3}, Values(rf), "RF1,RF2,RF3", "force", + "nodal")); + add(Data(root + "/nodal/residual", {n, 3}, Values(r), "R1,R2,R3", "force", + "nodal")); + D stress, strain; + Check(points.size() == 8 * e, "Solid point count is invalid."); + for (std::size_t i = 0; i < points.size(); ++i) { + Check(points[i].element == i / 8 && + points[i].integration_point == static_cast(i % 8 + 1), + "Solid point identity is invalid."); + stress.insert(stress.end(), points[i].cauchy_stress.begin(), + points[i].cauchy_stress.end()); + strain.insert(strain.end(), points[i].strain.begin(), + points[i].strain.end()); + } + auto s = + Data(root + "/solid/stress", {e, 8, 6}, std::move(stress), + "S11,S22,S33,S12,S13,S23", "force/length^2", "integration-point"); + s.strings["stress_measure"] = "cauchy"; + add(std::move(s)); + auto eps = + Data(root + "/solid/strain", {e, 8, 6}, std::move(strain), + "E11,E22,E33,E12,E13,E23", "dimensionless", "integration-point"); + eps.strings["strain_measure"] = domain.Steps()[0].UsesNonlinearGeometry() + ? "green-lagrange" + : "infinitesimal"; + eps.strings["shear_convention"] = "tensor"; + add(std::move(eps)); + auto en = Data(root + "/global/strain_energy", {}, D{energy}, "scalar", + "force*length", "global"); + en.strings["definition"] = domain.Steps()[0].UsesNonlinearGeometry() + ? "sum(W*detJ0)" + : "0.5*u^T*K0*u"; + add(std::move(en)); +} +void Model(SolidPlan& plan, const Domain& domain, const DofManager& dofs) { + SolidTable nodes{"/model/nodes", + domain.Nodes().size(), + {{"internal_node_id", 1, U{}}, + {"instance_name", 1, S{}}, + {"source_label", 1, S{}}, + {"coordinates", 3, D{}}}}; + for (std::size_t i = 0; i < domain.Nodes().size(); ++i) { + const auto& node = domain.Nodes()[i]; + std::get(nodes.columns[0].values).push_back(i); + std::get(nodes.columns[1].values) + .push_back(node.source_id.instance_name); + std::get(nodes.columns[2].values) + .push_back(node.source_id.source_label_text); + auto& coords = std::get(nodes.columns[3].values); + coords.insert(coords.end(), node.coordinates.begin(), + node.coordinates.end()); + } + plan.tables.push_back(std::move(nodes)); + SolidTable elements{"/model/elements", + domain.SolidElements().Size(), + {{"internal_element_id", 1, U{}}, + {"instance_name", 1, S{}}, + {"source_label", 1, S{}}, + {"source_element_type", 1, S{}}, + {"internal_formulation", 1, S{}}, + {"node_internal_ids", 8, U{}}, + {"solid_section_internal_id", 1, U{}}, + {"material_internal_id", 1, U{}}}}; + for (std::size_t i = 0; i < domain.SolidElements().Size(); ++i) { + const auto& element = domain.SolidElements()[i]; + std::get(elements.columns[0].values).push_back(i); + std::get(elements.columns[1].values) + .push_back(element.SourceId().instance_name); + std::get(elements.columns[2].values) + .push_back(element.SourceId().source_label_text); + std::get(elements.columns[3].values) + .push_back(std::string(element.SourceElementType())); + std::get(elements.columns[4].values).push_back("FESA-Hex8"); + auto& connectivity = std::get(elements.columns[5].values); + connectivity.insert(connectivity.end(), element.node_indices.begin(), + element.node_indices.end()); + std::get(elements.columns[6].values).push_back(element.section_index); + std::get(elements.columns[7].values).push_back(element.MaterialIndex()); + } + plan.tables.push_back(std::move(elements)); + SolidTable materials{"/model/solid/materials", + domain.Materials().Size(), + {{"internal_material_id", 1, U{}}, + {"name", 1, S{}}, + {"youngs_modulus", 1, D{}}, + {"poisson_ratio", 1, D{}}, + {"density_present", 1, B{}}, + {"density", 1, D{}}}}; + for (std::size_t i = 0; i < domain.Materials().Size(); ++i) { + const auto* material = + dynamic_cast(&domain.Materials()[i]); + Check(material != nullptr, "Solid material is not isotropic elastic."); + bool present = false; + double rho = 0; + for (const auto& density : domain.Densities()) + if (density.material_index == i) { + present = true; + rho = density.rho; + } + std::get(materials.columns[0].values).push_back(i); + std::get(materials.columns[1].values).push_back(material->Name()); + std::get(materials.columns[2].values) + .push_back(material->YoungsModulus()); + std::get(materials.columns[3].values) + .push_back(material->PoissonsRatio()); + std::get(materials.columns[4].values).push_back(present ? 1 : 0); + std::get(materials.columns[5].values).push_back(rho); + } + plan.tables.push_back(std::move(materials)); + SolidTable sections{"/model/solid/sections", + domain.SolidSections().Size(), + {{"internal_section_id", 1, U{}}, + {"source_file", 1, S{}}, + {"source_line", 1, U{}}, + {"source_elset", 1, S{}}, + {"material_internal_id", 1, U{}}}}; + for (std::size_t i = 0; i < domain.SolidSections().Size(); ++i) { + const auto& section = domain.SolidSections()[i]; + std::get(sections.columns[0].values).push_back(i); + std::get(sections.columns[1].values) + .push_back(section.Location().file.generic_u8string()); + std::get(sections.columns[2].values).push_back(section.Location().line); + std::get(sections.columns[3].values).push_back(section.Name()); + std::get(sections.columns[4].values).push_back(section.MaterialIndex()); + } + plan.tables.push_back(std::move(sections)); + SolidTable points{"/model/solid_integration_points", + 8 * domain.Elements().Size(), + {{"internal_element_id", 1, U{}}, + {"point_id", 1, B{}}, + {"natural_coordinates", 3, D{}}, + {"reference_coordinates", 3, D{}}}}; + ElementFactory factory; + for (std::size_t i = 0; i < domain.Elements().Size(); ++i) { + auto element = factory.Create(domain.Elements()[i], domain); + Check(element.HasValue(), "Solid factory reference failed."); + const auto* hex = dynamic_cast(element.Value().get()); + Check(hex != nullptr, "Solid kernel identity failed."); + auto response = hex->EvaluateLinear(Vector{24}); + Check(response.HasValue(), "Solid reference recovery failed."); + for (const auto& point : response.Value().points) { + std::get(points.columns[0].values).push_back(i); + std::get(points.columns[1].values) + .push_back(static_cast(point.integration_point)); + auto& natural = std::get(points.columns[2].values); + natural.insert(natural.end(), point.natural_coordinates.begin(), + point.natural_coordinates.end()); + auto& position = std::get(points.columns[3].values); + position.insert(position.end(), point.reference_coordinates.begin(), + point.reference_coordinates.end()); + } + } + plan.tables.push_back(std::move(points)); + plan.mask.assign(dofs.FullDofCount(), 0); + for (const auto index : dofs.ConstrainedDofs()) plan.mask[index] = 1; +} +void Diagnostics(SolidPlan& plan, std::vector diagnostics) { + SortDiagnostics(diagnostics); + SolidTable table{"/diagnostics", + diagnostics.size(), + {{"severity", 1, S{}}, + {"code", 1, S{}}, + {"file", 1, S{}}, + {"line", 1, U{}}, + {"keyword", 1, S{}}, + {"entity_identity", 1, S{}}, + {"message", 1, S{}}}}; + for (const auto& d : diagnostics) { + std::get(table.columns[0].values) + .push_back(d.severity == Severity::kWarning ? "warning" : "error"); + std::get(table.columns[1].values).push_back(d.code); + std::get(table.columns[2].values) + .push_back(d.location.file.generic_u8string()); + std::get(table.columns[3].values).push_back(d.location.line); + std::get(table.columns[4].values).push_back(d.keyword); + std::get(table.columns[5].values).push_back(d.entity_identity); + std::get(table.columns[6].values).push_back(d.message); + } + plan.tables.push_back(std::move(table)); +} +} // namespace + +Hdf5Handle SolidColumnType(const SolidColumn& column, bool native) { + Hdf5Handle result; + if (std::holds_alternative(column.values)) + result = MakeUtf8StringType(); + else { + hid_t type = std::holds_alternative(column.values) + ? (native ? H5T_NATIVE_UINT64 : H5T_STD_U64LE) + : std::holds_alternative(column.values) + ? (native ? H5T_NATIVE_UINT8 : H5T_STD_U8LE) + : (native ? H5T_NATIVE_DOUBLE : H5T_IEEE_F64LE); + result = Hdf5Handle{ + RequireHdf5Id(H5Tcopy(type), "Solid column type failed."), H5Tclose}; + } + if (column.width > 1) { + const hsize_t width = column.width; + result = Hdf5Handle{RequireHdf5Id(H5Tarray_create2(result.Get(), 1, &width), + "Solid array type failed."), + H5Tclose}; + } + return result; +} + +Result BuildSolidPlan(const Domain& domain, SolidStateView view, + const std::vector& diagnostics) { + try { + if (std::holds_alternative< + std::reference_wrapper>(view)) + Check(math_internal::ArithmeticEnvironmentSupported(), + "Solid RN64 serialization requires supported arithmetic " + "environment."); + Check(domain.IsSolid() && domain.Steps().Size() == 1, + "Solid output requires exactly one Solid step."); + auto model = AnalysisModel::Create(domain); + Check(model.HasValue(), "Solid output model failed."); + auto dofs = DofManager::Create(model.Value()); + Check(dofs.HasValue(), "Solid output DOF partition failed."); + SolidPlan plan; + const bool dynamic = std::holds_alternative< + std::reference_wrapper>(view); + const bool frequency = std::holds_alternative< + std::reference_wrapper>(view); + const bool nonlinear = domain.Steps()[0].UsesNonlinearGeometry(); + Check(dynamic == (domain.Dynamic() != nullptr) && + frequency == domain.Steps()[0].IsFrequency(), + "Solid output procedure/state mismatch."); + plan.metadata = { + {"feature_id", "hex8-structural-analysis"}, + {"analysis_procedure", dynamic ? (nonlinear ? "hht-nonlinear-dynamics" + : "hht-linear-dynamics") + : frequency ? "linear-frequency" + : nonlinear ? "nonlinear-static-full-newton" + : "linear-static"}, + {"element_formulation", "FESA-Hex8"}, + {"integration_rule", "gauss-2x2x2"}, + {"mass_formulation", "consistent-reference"}, + {"source_input_identity", SourceInputIdentity(domain)}, + {"source_step_name", domain.Steps()[0].Name()}, + {"unit_system_label", "user-consistent-unspecified"}, + {"coordinate_convention", + "global-cartesian; reference-strain; current-cauchy-stress"}, + {"component_order", "node-major; UX,UY,UZ"}, + {"solver_version", std::string(SolverVersion())}}; + Model(plan, domain, dofs.Value()); + Diagnostics(plan, diagnostics); + if (const auto* ref = + std::get_if>( + &view)) { + const auto& state = ref->get(); + auto validator = SolidStaticState::Create(domain, dofs.Value()); + Check(validator.HasValue(), "Static output validator creation failed."); + Check(!state.AcceptedFrames().empty(), "Static output history is empty."); + SolidTable accepted{"/steps/Step-1/accepted_frames", + state.AcceptedFrames().size(), + {{"frame_index", 1, U{}}, + {"increment_index", 1, U{}}, + {"load_factor", 1, D{}}, + {"iteration_count", 1, U{}}, + {"is_factor_one", 1, B{}}}}; + for (const auto& frame : state.AcceptedFrames()) { + Check(validator.Value().CommitAcceptedFrame(frame).IsOk(), + "Static accepted frame is invalid."); + FrameData(plan, domain, frame.frame_index, frame.increment_index, + frame.load_factor, false, frame.displacement, frame.reaction, + frame.residual, frame.solid_results, + frame.physical_strain_energy); + std::get(accepted.columns[0].values).push_back(frame.frame_index); + std::get(accepted.columns[1].values) + .push_back(frame.increment_index); + std::get(accepted.columns[2].values).push_back(frame.load_factor); + std::get(accepted.columns[3].values) + .push_back(frame.iteration_count); + std::get(accepted.columns[4].values) + .push_back(frame.is_factor_one ? 1 : 0); + } + Check(state.AcceptedFrames().back().is_factor_one, + "Static output is incomplete."); + if (nonlinear) plan.tables.push_back(std::move(accepted)); + } else if (const auto* dynamic_ref = std::get_if< + std::reference_wrapper>( + &view)) { + Check(math_internal::ArithmeticEnvironmentSupported(), + "Solid RN64 serialization requires supported arithmetic " + "environment."); + const auto& state = dynamic_ref->get(); + auto validator = SolidDynamicAnalysisState::Create(domain, dofs.Value()); + Check(validator.HasValue(), "Dynamic output validator creation failed."); + Check(state.AcceptedFrames().size() == + domain.Dynamic()->increment_count + 1, + "Dynamic output history is incomplete."); + SolidTable accepted{"/steps/Step-1/accepted_frames", + state.AcceptedFrames().size(), + {{"frame_index", 1, U{}}, + {"increment_index", 1, U{}}, + {"time", 1, D{}}, + {"iteration_count", 1, U{}}, + {"is_initial", 1, B{}}}}; + for (const auto& frame : state.AcceptedFrames()) { + Check(validator.Value().CommitAcceptedFrame(frame).IsOk(), + "Full-pair dynamic accepted frame is invalid."); + FrameData(plan, domain, frame.frame_index, frame.increment_index, + frame.time, true, frame.displacement, frame.reaction, + frame.residual, frame.solid_results, frame.strain_energy); + const std::string root = + "/steps/Step-1/frames/" + std::to_string(frame.frame_index); + const auto add = [&](SolidData data) { + data.integers = {{"frame_index", frame.frame_index}, + {"increment_index", frame.increment_index}}; + data.doubles = {{"time", frame.time}}; + plan.datasets.push_back(std::move(data)); + }; + add(Data(root + "/nodal/velocity", {domain.Nodes().size(), 3}, + Values(frame.velocity), "V1,V2,V3", "length/time", "nodal")); + auto a = Data(root + "/nodal/acceleration", {domain.Nodes().size(), 3}, + Values(frame.acceleration), "A1,A2,A3", "length/time^2", + "nodal"); + a.strings["acceleration_kind"] = "physical"; + add(std::move(a)); + a = Data(root + "/nodal/algorithmic_acceleration", + {domain.Nodes().size(), 3}, + Values(frame.algorithmic_acceleration), "A1,A2,A3", + "length/time^2", "nodal"); + a.strings["acceleration_kind"] = "algorithmic"; + add(std::move(a)); + for (const auto& quantity : + std::vector>{ + {"kinetic_energy", frame.kinetic_energy, "0.5*v^T*M*v"}, + {"external_work", frame.external_work, + "sum(0.5*(F_old+F_new)*delta_u)"}, + {"energy_balance_error", frame.energy_balance_error, + "T+U-T0-U0-Wext"}}) { + auto data = Data(root + "/global/" + std::get<0>(quantity), {}, + D{std::get<1>(quantity)}, "scalar", "force*length", + "global"); + data.strings["definition"] = std::get<2>(quantity); + add(std::move(data)); + } + add(Data(root + "/global/linear_momentum", {3}, + D(frame.linear_momentum.begin(), frame.linear_momentum.end()), + "P1,P2,P3", "mass*length/time", "global")); + auto angular = Data( + root + "/global/angular_momentum", {3}, + D(frame.angular_momentum.begin(), frame.angular_momentum.end()), + "L1,L2,L3", "mass*length^2/time", "global"); + angular.strings["origin"] = "global0"; + add(std::move(angular)); + std::get(accepted.columns[0].values).push_back(frame.frame_index); + std::get(accepted.columns[1].values) + .push_back(frame.increment_index); + std::get(accepted.columns[2].values).push_back(frame.time); + std::get(accepted.columns[3].values) + .push_back(frame.iteration_count); + std::get(accepted.columns[4].values) + .push_back(frame.is_initial ? 1 : 0); + } + plan.tables.push_back(std::move(accepted)); + } else { + const auto& state = + std::get>(view) + .get(); + Check(domain.Frequency()->requested_mode_count > 0, + "Modal request is invalid."); + const std::size_t count = + static_cast(domain.Frequency()->requested_mode_count); + Check(state.full_dof_count == dofs.Value().FullDofCount() && + state.free_dof_count == dofs.Value().FreeDofs().size() && + state.constraint_mask == plan.mask, + "Modal DOF inventory is invalid."); + Check(state.prescribed_displacement.size() == plan.mask.size(), + "Modal prescribed shape is invalid."); + for (const auto value : state.prescribed_displacement) + Check(value == 0, "Modal prescribed value is invalid."); + Check(state.modes.size() == count && state.vectors.Rows() == count && + state.vectors.Columns() == plan.mask.size(), + "Modal output shape is invalid."); + Check(std::isfinite(state.classification_scale) && + state.classification_scale > 0 && + std::isfinite(state.mass_lower_bound) && + state.mass_lower_bound > 0, + "Modal certification scale is invalid."); + Check(std::isfinite(state.eta_gram) && state.eta_gram >= 0 && + state.eta_gram <= 1e-10 && std::isfinite(state.max_eta_orth) && + state.max_eta_orth >= 0 && state.max_eta_orth <= 1e-10, + "Modal Gram evidence is invalid."); + Check(state.eigenpair_checks.size() == state.backend.eigenvalues.size() && + state.backend.column_ids.size() == + state.eigenpair_checks.size() && + state.eigenpair_checks.size() >= count && + state.eigenpair_checks.size() <= state.free_dof_count, + "Modal backend verification inventory is invalid."); + std::size_t expected_cluster = 1; + SolidTable modes{"/steps/Step-1/eigen/modes", + count, + {{"mode_id", 1, U{}}, + {"eigenvalue", 1, D{}}, + {"classification", 1, B{}}, + {"cluster_id", 1, U{}}}}; + SolidTable frequencies{"/steps/Step-1/eigen/frequencies", + 0, + {{"mode_id", 1, U{}}, {"frequency", 1, D{}}}}; + D vectors; + for (std::size_t i = 0; i < count; ++i) { + const auto& mode = state.modes[i]; + Check(mode.mode_id == i + 1 && mode.cluster_id > 0 && + std::isfinite(mode.eigenvalue) && + mode.eigenvalue / state.classification_scale >= -1e-9 && + (i == 0 || mode.eigenvalue >= state.modes[i - 1].eigenvalue), + "Modal identity/order is invalid."); + if (i > 0 && (mode.eigenvalue - state.modes[i - 1].eigenvalue) / + state.classification_scale > + 1e-8) + ++expected_cluster; + Check(mode.cluster_id == expected_cluster && + mode.eigenvalue == state.eigenpair_checks[i].eigenvalue, + "Modal cluster or selected-pair identity is invalid."); + const bool near = + std::abs(mode.eigenvalue) / state.classification_scale <= 1e-9; + Check(mode.classification == (near ? 0 : 1), + "Modal classification is invalid."); + Check(mode.frequency.has_value() == (mode.eigenvalue >= 0), + "Modal frequency presence is invalid."); + if (mode.frequency) { + Check(*mode.frequency == + std::sqrt(mode.eigenvalue) / (2 * std::acos(-1.0)), + "Modal frequency is invalid."); + ++frequencies.rows; + std::get(frequencies.columns[0].values).push_back(mode.mode_id); + std::get(frequencies.columns[1].values).push_back(*mode.frequency); + } + std::get(modes.columns[0].values).push_back(mode.mode_id); + std::get(modes.columns[1].values).push_back(mode.eigenvalue); + std::get(modes.columns[2].values).push_back(mode.classification); + std::get(modes.columns[3].values).push_back(mode.cluster_id); + for (std::size_t j = 0; j < plan.mask.size(); ++j) { + Check(std::isfinite(state.vectors(i, j)) && + (!plan.mask[j] || state.vectors(i, j) == 0), + "Modal vector is invalid."); + vectors.push_back(state.vectors(i, j)); + } + } + modes.strings = {{"component_names", "eigenvalue"}, + {"unit_dimensions", "1/time^2"}, + {"location", "global"}, + {"coordinate_system", "global-cartesian"}}; + frequencies.strings = {{"component_names", "frequency"}, + {"unit_dimensions", "1/time"}, + {"location", "global"}, + {"coordinate_system", "global-cartesian"}}; + plan.tables.push_back(std::move(modes)); + plan.tables.push_back(std::move(frequencies)); + auto data = + Data("/steps/Step-1/eigen/vectors", {count, domain.Nodes().size(), 3}, + std::move(vectors), "UX,UY,UZ", "mass-normalized", "nodal"); + data.strings["normalization"] = "phi^T*M*phi=1"; + plan.datasets.push_back(std::move(data)); + SolidTable pairs{"/steps/Step-1/eigen/verification/pairs", + state.eigenpair_checks.size(), + {{"sorted_rank", 1, U{}}, + {"solver_column", 1, U{}}, + {"eigenvalue", 1, D{}}, + {"modal_mass", 1, D{}}, + {"rayleigh_quotient", 1, D{}}, + {"eta_r", 1, D{}}, + {"eta_norm", 1, D{}}, + {"eta_energy", 1, D{}}}}; + + std::set solver_columns; + for (std::size_t i = 0; i < state.eigenpair_checks.size(); ++i) { + Check(state.backend.column_ids[i] < state.eigenpair_checks.size() && + solver_columns.insert(state.backend.column_ids[i]).second, + "Modal backend column inventory is invalid."); + const auto& pair = state.eigenpair_checks[i]; + const double normalized = pair.eigenvalue / state.classification_scale; + Check(std::isfinite(pair.eigenvalue) && std::isfinite(normalized), + "Backend eigenvalue or normalized eigenvalue is nonfinite."); + Check(normalized >= -1e-9, + "Significant negative backend eigenvalue cannot be published."); + Check(i == 0 || + pair.eigenvalue >= state.eigenpair_checks[i - 1].eigenvalue, + "Backend eigenvalues are not sorted."); + + Check(pair.sorted_rank == i + 1 && + pair.solver_column == state.backend.column_ids[i] && + pair.eigenvalue == state.backend.eigenvalues[i], + "Modal pair identity is invalid."); + Check(std::isfinite(pair.modal_mass) && pair.modal_mass > 0 && + std::abs(pair.modal_mass - 1) <= 1e-10, + "Modal mass evidence is invalid."); + for (const double eta : {pair.eta_r, pair.eta_norm, pair.eta_energy}) + Check(std::isfinite(eta) && eta >= 0 && eta <= 1e-10, + "Modal pair gate is invalid."); + std::get(pairs.columns[0].values).push_back(pair.sorted_rank); + std::get(pairs.columns[1].values).push_back(pair.solver_column); + const std::array values{{pair.eigenvalue, pair.modal_mass, + pair.rayleigh_quotient, pair.eta_r, + pair.eta_norm, pair.eta_energy}}; + for (std::size_t j = 0; j < 6; ++j) + std::get(pairs.columns[j + 2].values).push_back(values[j]); + } + plan.tables.push_back(std::move(pairs)); + plan.datasets.push_back(Data("/steps/Step-1/eigen/verification/eta_gram", + {}, D{state.eta_gram}, "scalar", + "dimensionless", "global")); + plan.datasets.push_back( + Data("/steps/Step-1/eigen/verification/max_eta_orth", {}, + D{state.max_eta_orth}, "scalar", "dimensionless", "global")); + } + for (const auto& table : plan.tables) + for (const auto& column : table.columns) + std::visit( + [&](const auto& values) { + Check(values.size() == table.rows * column.width, + "Solid column shape is invalid."); + using T = typename std::decay_t::value_type; + for (const auto& value : values) { + if constexpr (std::is_same_v) + Check(std::isfinite(value), + "Solid numeric column is nonfinite."); + if constexpr (std::is_same_v) + Check(IsValidUtf8(value) && + value.find('\0') == std::string::npos, + "Solid string column is invalid UTF8."); + } + }, + column.values); + for (const auto& data : plan.datasets) { + std::size_t size = 1; + for (const auto width : data.dimensions) { + Check(SizeProductFits(size, width, size), + "Solid dataset shape overflow."); + } + Check(size == data.values.size(), "Solid dataset shape is invalid."); + for (const auto value : data.values) + Check(std::isfinite(value), "Solid dataset value is nonfinite."); + } + for (const auto& item : plan.metadata) + Check(IsValidUtf8(item.second) && + item.second.find('\0') == std::string::npos, + "Solid metadata UTF8 is invalid."); + return Result::Success(std::move(plan)); + } catch (const std::exception& error) { + return Result::Failure( + OutputFailure("invalid-solid-output", error.what())); + } +} + +Status WriteSolidFile(const std::filesystem::path& output, const Domain& domain, + SolidStateView view, + const std::vector& diagnostics, + const SolidFileHooks& hooks) { + const auto plan = BuildSolidPlan(domain, view, diagnostics); + if (!plan.HasValue()) return plan.GetStatus(); + try { + Hdf5ErrorSilencer silence; + const auto candidate = MakeCandidatePath(output); + CandidateFileGuard guard(candidate); + Hdf5Handle file{ + RequireHdf5Id(H5Fcreate(candidate.string().c_str(), H5F_ACC_EXCL, + H5P_DEFAULT, H5P_DEFAULT), + "Solid candidate create failed."), + hooks.candidate_close}; + if (hooks.before_write) hooks.before_write(file.Get()); + if (std::holds_alternative< + std::reference_wrapper>(view)) + Check(math_internal::ArithmeticEnvironmentSupported(), + "Solid RN64 serialization arithmetic environment changed."); + { + auto group = CreateGroup(file.Get(), "/metadata"); + WriteUint64Attribute(group.Get(), "schema_version", 1); + for (const auto& item : plan.Value().metadata) + WriteStringAttribute(group.Get(), item.first.c_str(), item.second); + } + (void)CreateGroup(file.Get(), "/model"); + WriteNodes(file.Get(), domain); + for (const auto& table : plan.Value().tables) + if (table.path != "/model/nodes" && table.path != "/diagnostics") + Table(file.Get(), table); + WriteDiagnostics(file.Get(), diagnostics); + { + const std::vector dimensions{domain.Nodes().size(), 3}; + auto space = CreateDatasetSpace(dimensions); + Hdf5Handle dataset{ + RequireHdf5Id(H5Dcreate2(file.Get(), "/model/nodal_constraint_mask", + H5T_STD_U8LE, space.Get(), H5P_DEFAULT, + H5P_DEFAULT, H5P_DEFAULT), + "Solid mask create failed."), + H5Dclose}; + RequireHdf5(H5Dwrite(dataset.Get(), H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL, + H5P_DEFAULT, plan.Value().mask.data()), + "Solid mask write failed."); + WriteStringAttribute(dataset.Get(), "component_names", "UX,UY,UZ"); + } + for (const auto& data : plan.Value().datasets) { + Parent(file.Get(), data.path); + auto dataset = WriteDoubleValues(file.Get(), data.path, data.dimensions, + data.values.data(), data.values.size()); + for (const auto& item : data.strings) + WriteStringAttribute(dataset.Get(), item.first.c_str(), item.second); + for (const auto& item : data.integers) + WriteUint64Attribute(dataset.Get(), item.first.c_str(), item.second); + for (const auto& item : data.doubles) + WriteDoubleAttribute(dataset.Get(), item.first.c_str(), item.second); + } + RequireHdf5(H5Fflush(file.Get(), H5F_SCOPE_GLOBAL), + "Solid candidate flush failed."); + RequireHdf5(file.CloseChecked(), "Solid candidate close failed."); + if (hooks.before_self_check) hooks.before_self_check(candidate); + const Status check = + SelfCheckSolidFile(candidate, plan.Value(), hooks.reader_close); + if (!check.IsOk()) return check; + if (hooks.before_finalize) hooks.before_finalize(candidate); + if (std::holds_alternative< + std::reference_wrapper>(view)) + Check(math_internal::ArithmeticEnvironmentSupported(), + "Solid RN64 publication arithmetic environment changed."); + const Status finalized = + AtomicFinalizeValidatedCandidate(candidate, output, check); + if (finalized.IsOk()) guard.Release(); + return finalized; + } catch (const std::exception& error) { + return OutputFailure("solid-output-failure", error.what()); + } +} +} // namespace fesa::hdf5_internal + +namespace fesa { +Status SolidHdf5ResultsWriter::Write( + const std::filesystem::path& path, const Domain& domain, + const SolidStaticState& state, const std::vector& diagnostics) { + return hdf5_internal::WriteSolidFile(path, domain, std::cref(state), + diagnostics); +} +Status SolidHdf5ResultsWriter::Write( + const std::filesystem::path& path, const Domain& domain, + const SolidFrequencyState& state, + const std::vector& diagnostics) { + return hdf5_internal::WriteSolidFile(path, domain, std::cref(state), + diagnostics); +} +Status SolidHdf5ResultsWriter::Write( + const std::filesystem::path& path, const Domain& domain, + const SolidDynamicAnalysisState& state, + const std::vector& diagnostics) { + return hdf5_internal::WriteSolidFile(path, domain, std::cref(state), + diagnostics); +} +} // namespace fesa diff --git a/src/fesa/io/hdf5/solid_hdf5_self_check.cpp b/src/fesa/io/hdf5/solid_hdf5_self_check.cpp new file mode 100644 index 0000000..0128150 --- /dev/null +++ b/src/fesa/io/hdf5/solid_hdf5_self_check.cpp @@ -0,0 +1,217 @@ +#include "io/hdf5/solid_hdf5_self_check.h" + +#include +#include + +#include "io/hdf5/hdf5_primitives.h" + +namespace fesa::hdf5_internal { +namespace { +void Check(bool condition, const char* message) { + if (!condition) throw Hdf5Failure(message); +} +Hdf5Handle Dataset(hid_t file, const std::string& path, + const std::vector& dimensions) { + Hdf5Handle dataset{RequireHdf5Id(H5Dopen2(file, path.c_str(), H5P_DEFAULT), + "Solid required dataset is missing."), + H5Dclose}; + Hdf5Handle space{ + RequireHdf5Id(H5Dget_space(dataset.Get()), "Solid dataset space failed."), + H5Sclose}; + Check(H5Sget_simple_extent_ndims(space.Get()) == + static_cast(dimensions.size()), + "Solid dataset rank mismatch."); + std::vector actual(dimensions.size()); + RequireHdf5(H5Sget_simple_extent_dims(space.Get(), actual.data(), nullptr), + "Solid dimensions read failed."); + Check(actual == dimensions, "Solid dataset shape mismatch."); + Hdf5Handle properties{RequireHdf5Id(H5Dget_create_plist(dataset.Get()), + "Solid dataset properties failed."), + H5Pclose}; + Check(H5Pget_layout(properties.Get()) == H5D_CONTIGUOUS && + H5Pget_nfilters(properties.Get()) == 0, + "Solid dataset storage mismatch."); + return dataset; +} +void String(hid_t object, const std::string& name, + const std::string& expected) { + Hdf5Handle attr{RequireHdf5Id(H5Aopen(object, name.c_str(), H5P_DEFAULT), + "Solid string attribute missing."), + H5Aclose}; + Hdf5Handle type{ + RequireHdf5Id(H5Aget_type(attr.Get()), "Solid attribute type failed."), + H5Tclose}; + Hdf5Handle space{ + RequireHdf5Id(H5Aget_space(attr.Get()), "Solid attribute space failed."), + H5Sclose}; + Check(H5Tget_class(type.Get()) == H5T_STRING && + H5Tis_variable_str(type.Get()) > 0 && + H5Tget_cset(type.Get()) == H5T_CSET_UTF8 && + H5Sget_simple_extent_type(space.Get()) == H5S_SCALAR, + "Solid string attribute type mismatch."); + char* value = nullptr; + RequireHdf5(H5Aread(attr.Get(), type.Get(), &value), + "Solid string attribute read failed."); + const std::string actual = value ? value : ""; + RequireHdf5(H5free_memory(value), "Solid string release failed."); + Check(IsValidUtf8(actual) && actual == expected, + "Solid string attribute mismatch."); +} +template +void Number(hid_t object, const std::string& name, T expected, hid_t ft, + hid_t mt) { + Hdf5Handle attr{RequireHdf5Id(H5Aopen(object, name.c_str(), H5P_DEFAULT), + "Solid numeric attribute missing."), + H5Aclose}; + Hdf5Handle type{ + RequireHdf5Id(H5Aget_type(attr.Get()), "Solid attribute type failed."), + H5Tclose}; + Hdf5Handle space{ + RequireHdf5Id(H5Aget_space(attr.Get()), "Solid attribute space failed."), + H5Sclose}; + Check(H5Tequal(type.Get(), ft) > 0 && + H5Sget_simple_extent_type(space.Get()) == H5S_SCALAR, + "Solid numeric attribute type mismatch."); + T actual{}; + RequireHdf5(H5Aread(attr.Get(), mt, &actual), + "Solid numeric attribute read failed."); + Check(actual == expected, "Solid numeric attribute mismatch."); +} +void Column(hid_t dataset, hid_t compound, std::size_t index, + const SolidColumn& column, std::size_t rows) { + char* name = H5Tget_member_name(compound, static_cast(index)); + Check(name != nullptr, "Solid member name missing."); + const std::string actual_name = name; + RequireHdf5(H5free_memory(name), "Solid member name release failed."); + Check(actual_name == column.name, + "Solid compound member order/name mismatch."); + Hdf5Handle stored{ + RequireHdf5Id(H5Tget_member_type(compound, static_cast(index)), + "Solid member type missing."), + H5Tclose}; + const auto expected = SolidColumnType(column, false); + Check(H5Tequal(stored.Get(), expected.Get()) > 0, + "Solid compound member type mismatch."); + const auto native = SolidColumnType(column, true); + Hdf5Handle memory{ + RequireHdf5Id(H5Tcreate(H5T_COMPOUND, H5Tget_size(native.Get())), + "Solid column memory type failed."), + H5Tclose}; + RequireHdf5(H5Tinsert(memory.Get(), column.name.c_str(), 0, native.Get()), + "Solid column memory member failed."); + std::visit( + [&](const auto& values) { + using T = typename std::decay_t::value_type; + if constexpr (std::is_same_v) { + std::vector actual(rows, nullptr); + if (rows) + RequireHdf5(H5Dread(dataset, memory.Get(), H5S_ALL, H5S_ALL, + H5P_DEFAULT, actual.data()), + "Solid string column read failed."); + bool same = true; + for (std::size_t i = 0; i < rows; ++i) { + const std::string value = actual[i] ? actual[i] : ""; + same = same && IsValidUtf8(value) && value == values[i]; + } + Hdf5Handle space{RequireHdf5Id(H5Dget_space(dataset), + "Solid reclaim space failed."), + H5Sclose}; + if (rows) + RequireHdf5(H5Dvlen_reclaim(memory.Get(), space.Get(), H5P_DEFAULT, + actual.data()), + "Solid string column release failed."); + Check(same, "Solid string column identity mismatch."); + } else { + std::vector actual(values.size()); + if (rows) + RequireHdf5(H5Dread(dataset, memory.Get(), H5S_ALL, H5S_ALL, + H5P_DEFAULT, actual.data()), + "Solid numeric column read failed."); + Check(actual == values, "Solid numeric column identity mismatch."); + if constexpr (std::is_same_v) + for (const auto value : actual) + Check(std::isfinite(value), "Solid numeric column is nonfinite."); + } + }, + column.values); +} +} // namespace +Status SelfCheckSolidFile(const std::filesystem::path& path, + const SolidPlan& plan, + herr_t (*reader_close)(hid_t)) { + try { + Hdf5ErrorSilencer silence; + Hdf5Handle file{RequireHdf5Id(H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, + H5P_DEFAULT), + "Solid candidate readback open failed."), + reader_close}; + { + Hdf5Handle group{ + RequireHdf5Id(H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), + "Solid metadata missing."), + H5Gclose}; + Number(group.Get(), "schema_version", 1, H5T_STD_U64LE, + H5T_NATIVE_UINT64); + for (const auto& item : plan.metadata) + String(group.Get(), item.first, item.second); + } + for (const auto& table : plan.tables) { + auto dataset = Dataset(file.Get(), table.path, {table.rows}); + Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), + "Solid compound type failed."), + H5Tclose}; + Check(H5Tget_class(type.Get()) == H5T_COMPOUND && + H5Tget_nmembers(type.Get()) == + static_cast(table.columns.size()), + "Solid compound member inventory mismatch."); + for (std::size_t i = 0; i < table.columns.size(); ++i) + Column(dataset.Get(), type.Get(), i, table.columns[i], table.rows); + for (const auto& item : table.strings) + String(dataset.Get(), item.first, item.second); + } + { + auto dataset = Dataset(file.Get(), "/model/nodal_constraint_mask", + {plan.mask.size() / 3, 3}); + Hdf5Handle type{ + RequireHdf5Id(H5Dget_type(dataset.Get()), "Solid mask type failed."), + H5Tclose}; + Check(H5Tequal(type.Get(), H5T_STD_U8LE) > 0, + "Solid mask dtype mismatch."); + std::vector mask(plan.mask.size()); + RequireHdf5(H5Dread(dataset.Get(), H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL, + H5P_DEFAULT, mask.data()), + "Solid mask read failed."); + Check(mask == plan.mask, "Solid mask identity mismatch."); + String(dataset.Get(), "component_names", "UX,UY,UZ"); + } + for (const auto& data : plan.datasets) { + auto dataset = Dataset(file.Get(), data.path, data.dimensions); + Hdf5Handle type{RequireHdf5Id(H5Dget_type(dataset.Get()), + "Solid result dtype failed."), + H5Tclose}; + Check(H5Tequal(type.Get(), H5T_IEEE_F64LE) > 0, + "Solid result dtype mismatch."); + std::vector values(data.values.size()); + if (!values.empty()) + RequireHdf5(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), + "Solid result read failed."); + Check(values == data.values, "Solid result exact readback mismatch."); + for (const auto value : values) + Check(std::isfinite(value), "Solid result is nonfinite."); + for (const auto& item : data.strings) + String(dataset.Get(), item.first, item.second); + for (const auto& item : data.integers) + Number(dataset.Get(), item.first, item.second, H5T_STD_U64LE, + H5T_NATIVE_UINT64); + for (const auto& item : data.doubles) + Number(dataset.Get(), item.first, item.second, H5T_IEEE_F64LE, + H5T_NATIVE_DOUBLE); + } + RequireHdf5(file.CloseChecked(), "Solid readback close failed."); + return Status::Ok(); + } catch (const std::exception& error) { + return OutputFailure("solid-self-check-failure", error.what()); + } +} +} // namespace fesa::hdf5_internal diff --git a/src/fesa/io/hdf5/solid_hdf5_self_check.h b/src/fesa/io/hdf5/solid_hdf5_self_check.h new file mode 100644 index 0000000..c28e3bc --- /dev/null +++ b/src/fesa/io/hdf5/solid_hdf5_self_check.h @@ -0,0 +1,64 @@ +#ifndef FESA_IO_HDF5_SOLID_HDF5_SELF_CHECK_H_ +#define FESA_IO_HDF5_SOLID_HDF5_SELF_CHECK_H_ +#include + +#include +#include +#include + +#include "fesa/results/solid_results_writer.h" +#include "io/hdf5/hdf5_raii.h" +namespace fesa::hdf5_internal { +/// @brief Exact prescribed Solid compound columns; width>1 denotes a fixed +/// array. +struct SolidColumn { + std::string name; + std::size_t width{1}; + std::variant, std::vector, + std::vector, std::vector> + values; +}; +struct SolidTable { + std::string path; + std::size_t rows{0}; + std::vector columns; + std::map strings; +}; +struct SolidData { + std::string path; + std::vector dimensions; + std::vector values; + std::map strings; + std::map integers; + std::map doubles; +}; +/// @brief Expected RN64 serialization inventory, never a rounded-only physics +/// check. +struct SolidPlan { + std::map metadata; + std::vector mask; + std::vector tables; + std::vector datasets; +}; +/// @brief Private fault seams around real writes/readback/close/replace. +struct SolidFileHooks { + herr_t (*candidate_close)(hid_t) = H5Fclose; + herr_t (*reader_close)(hid_t) = H5Fclose; + std::function before_write; + std::function before_self_check; + std::function before_finalize; +}; +using SolidStateView = + std::variant, + std::reference_wrapper, + std::reference_wrapper>; +Result BuildSolidPlan(const Domain&, SolidStateView, + const std::vector&); +Status WriteSolidFile(const std::filesystem::path&, const Domain&, + SolidStateView, const std::vector&, + const SolidFileHooks& = {}); +Status SelfCheckSolidFile(const std::filesystem::path&, const SolidPlan&, + herr_t (*reader_close)(hid_t) = H5Fclose); +Hdf5Handle SolidColumnType(const SolidColumn&, bool native); +} // namespace fesa::hdf5_internal +#endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f1a270e..cab0209 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -71,6 +71,7 @@ add_executable( unit/io/hdf5/nonlinear_hdf5_results_writer_test.cpp unit/io/hdf5/hdf5_writer_components_test.cpp unit/io/hdf5/thermal_hdf5_results_writer_test.cpp + unit/io/hdf5/solid_hdf5_results_writer_test.cpp unit/loads/load_test.cpp unit/loads/shell_surface_pressure_test.cpp unit/model/domain_test.cpp @@ -118,6 +119,7 @@ add_executable( integration/app/eigen_application_test.cpp integration/app/shell_surface_pressure_application_test.cpp integration/app/thermal_application_test.cpp + integration/app/solid_application_test.cpp ) add_dependencies(fesa_integration_tests fesa_cli) diff --git a/tests/integration/app/solid_application_test.cpp b/tests/integration/app/solid_application_test.cpp new file mode 100644 index 0000000..5163581 --- /dev/null +++ b/tests/integration/app/solid_application_test.cpp @@ -0,0 +1,174 @@ +#include +#include + +#include +#include +#include +#include +#include +#define NOMINMAX +#include + +#include "../../unit/io/hdf5/solid_output_test_fixture.h" +#include "fesa/app/fesa_application.h" +namespace { + +std::map> NumericSnapshot( + const std::filesystem::path& path, unsigned p) { + std::vector paths{"/model/nodal_constraint_mask"}; + if (p == 2) + paths.insert( + paths.end(), + {"/steps/Step-1/eigen/modes", "/steps/Step-1/eigen/frequencies", + "/steps/Step-1/eigen/vectors", + "/steps/Step-1/eigen/verification/pairs", + "/steps/Step-1/eigen/verification/eta_gram", + "/steps/Step-1/eigen/verification/max_eta_orth"}); + else { + if (p != 0) paths.push_back("/steps/Step-1/accepted_frames"); + const unsigned count = p == 0 ? 1 : p == 1 ? 2 : 3; + for (unsigned i = 0; i < count; ++i) { + const auto root = "/steps/Step-1/frames/" + std::to_string(i); + for (const auto* field : + {"nodal/displacement", "nodal/reaction", "nodal/residual", + "solid/stress", "solid/strain", "global/strain_energy"}) + paths.push_back(root + "/" + field); + if (p >= 3) + for (const auto* field : + {"nodal/velocity", "nodal/acceleration", + "nodal/algorithmic_acceleration", "global/kinetic_energy", + "global/external_work", "global/energy_balance_error", + "global/linear_momentum", "global/angular_momentum"}) + paths.push_back(root + "/" + field); + } + } + const hid_t file = + H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + if (file < 0) throw std::runtime_error("CLI result file missing"); + std::map> result; + for (const auto& dataset_path : paths) { + const hid_t dataset = H5Dopen2(file, dataset_path.c_str(), H5P_DEFAULT); + if (dataset < 0) { + H5Fclose(file); + throw std::runtime_error("CLI mandatory dataset missing"); + } + const hid_t type = H5Dget_type(dataset), space = H5Dget_space(dataset); + const auto count = H5Sget_simple_extent_npoints(space); + std::vector bytes(static_cast(count) * + H5Tget_size(type)); + const herr_t read = count == 0 ? 0 + : H5Dread(dataset, type, H5S_ALL, H5S_ALL, + H5P_DEFAULT, bytes.data()); + H5Sclose(space); + H5Tclose(type); + H5Dclose(dataset); + if (read < 0) { + H5Fclose(file); + throw std::runtime_error("CLI mandatory dataset read failed"); + } + result.emplace(dataset_path, std::move(bytes)); + } + H5Fclose(file); + return result; +} +TEST(SolidApplication, AllFiveActualProceduresWriteOfficialSolidSchema) { + for (unsigned p = 0; p < 5; ++p) { + SCOPED_TRACE(p); + solid_output_test::Files files(p); + ASSERT_EQ(fesa::FesaApplication{}.Run({files.input.string()}), 0); + const hid_t f = + H5Fopen(files.output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + ASSERT_GE(f, 0); + const hid_t g = H5Gopen2(f, "/metadata", H5P_DEFAULT); + ASSERT_GE(g, 0); + const hid_t a = H5Aopen(g, "analysis_procedure", H5P_DEFAULT); + ASSERT_GE(a, 0); + const hid_t t = H5Aget_type(a); + char* value = nullptr; + ASSERT_GE(H5Aread(a, t, &value), 0); + EXPECT_STREQ(value, solid_output_test::Procedure(p)); + H5free_memory(value); + H5Tclose(t); + H5Aclose(a); + H5Gclose(g); + H5Fclose(f); + } +} + +TEST(SolidApplication, + NativeExecutableRunsAllFiveAndKeepsWarningsAndTypedLogs) { + std::vector module(32768); + const DWORD length = GetModuleFileNameW(nullptr, module.data(), + static_cast(module.size())); + ASSERT_GT(length, 0u); + ASSERT_LT(length, module.size()); + const auto executable = + std::filesystem::path(std::wstring(module.data(), length)) + .parent_path() + .parent_path() / + "src/fesa/fesa.exe"; + ASSERT_TRUE(std::filesystem::exists(executable)); + for (unsigned p = 0; p < 5; ++p) { + SCOPED_TRACE(p); + solid_output_test::Files files(p); + ASSERT_EQ(fesa::FesaApplication{}.Run({files.input.string()}), 0); + const auto expected = NumericSnapshot(files.output, p); + std::wstring command = + L"\"" + executable.wstring() + L"\" \"" + files.input.wstring() + L"\""; + STARTUPINFOW startup{}; + startup.cb = sizeof(startup); + PROCESS_INFORMATION process{}; + ASSERT_TRUE(CreateProcessW(executable.c_str(), command.data(), nullptr, + nullptr, FALSE, CREATE_NO_WINDOW, nullptr, + nullptr, &startup, &process)); + const DWORD wait = WaitForSingleObject(process.hProcess, 60000); + if (wait != WAIT_OBJECT_0) TerminateProcess(process.hProcess, 99); + DWORD exit = 99; + EXPECT_TRUE(GetExitCodeProcess(process.hProcess, &exit)); + CloseHandle(process.hThread); + CloseHandle(process.hProcess); + ASSERT_EQ(wait, WAIT_OBJECT_0); + ASSERT_EQ(exit, 0u); + EXPECT_GT(H5Fis_hdf5(files.output.string().c_str()), 0); + EXPECT_EQ(NumericSnapshot(files.output, p), expected); + std::ifstream stream(files.log, std::ios::binary); + std::string log{std::istreambuf_iterator(stream), {}}; + EXPECT_EQ(log.find("correction_rotation"), std::string::npos); + EXPECT_EQ(log.find("residual_moment"), std::string::npos); + if (p == 1) + EXPECT_NE(log.find("schema=solid-static-log-v1"), std::string::npos); + if (p >= 3) { + EXPECT_NE(log.find("schema=solid-hht-log-v1"), std::string::npos); + EXPECT_NE(log.find("h2_physical_da_over_L="), std::string::npos); + EXPECT_NE(log.find("h2_algorithmic_dA_over_L="), std::string::npos); + EXPECT_NE(log.find("normalized-high-low-physical-state;RN64-observables"), + std::string::npos); + } + const hid_t file = + H5Fopen(files.output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + ASSERT_GE(file, 0); + const hid_t diagnostics = H5Dopen2(file, "/diagnostics", H5P_DEFAULT); + ASSERT_GE(diagnostics, 0); + const hid_t space = H5Dget_space(diagnostics); + hsize_t rows = 0; + H5Sget_simple_extent_dims(space, &rows, nullptr); + EXPECT_GT(rows, 0u); + H5Sclose(space); + H5Dclose(diagnostics); + H5Fclose(file); + } +} +TEST(SolidApplication, SuccessfulWarningsRetainSourceLocationAndTarget) { + solid_output_test::Files files(0); + testing::internal::CaptureStderr(); + const int exit = fesa::FesaApplication{}.Run({files.input.string()}); + const auto diagnostics = testing::internal::GetCapturedStderr(); + ASSERT_EQ(exit, 0); + EXPECT_NE(diagnostics.find("code=inactive-solid-rotation-boundary"), + std::string::npos); + EXPECT_NE(diagnostics.find("entity_identity=Fixed"), std::string::npos); + EXPECT_NE(diagnostics.find("line=30"), std::string::npos); + EXPECT_NE(diagnostics.find(files.input.generic_u8string()), + std::string::npos); +} +} // namespace diff --git a/tests/unit/io/hdf5/solid_hdf5_results_writer_test.cpp b/tests/unit/io/hdf5/solid_hdf5_results_writer_test.cpp new file mode 100644 index 0000000..84235d5 --- /dev/null +++ b/tests/unit/io/hdf5/solid_hdf5_results_writer_test.cpp @@ -0,0 +1,871 @@ +#include "fesa/io/hdf5/solid_hdf5_results_writer.h" + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "fesa/analysis/linear_frequency_procedure.h" +#include "fesa/analysis/linear_static_procedure.h" +#include "fesa/analysis/newton_iteration_log.h" +#include "fesa/analysis/nonlinear_static_procedure.h" +#include "fesa/analysis/prepared_analysis.h" +#include "fesa/analysis/solid_frequency_procedure.h" +#include "fesa/app/output_paths.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/io/hdf5/eigen_hdf5_results_writer.h" +#include "fesa/io/hdf5/hdf5_results_writer.h" +#include "fesa/io/hdf5/nonlinear_hdf5_results_writer.h" +#include "fesa/solvers/linear/mkl_pardiso_solver.h" +#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" +#include "io/hdf5/hdf5_primitives.h" +#include "io/hdf5/solid_hdf5_self_check.h" +#include "solid_output_test_fixture.h" +namespace { +using namespace fesa; +using namespace fesa::hdf5_internal; +struct Actual { + PreparedAnalysisInput input; + PreparedAnalysisOutcome outcome; +}; +Actual Execute(const solid_output_test::Files& files) { + auto input = PrepareAnalysisInput(files.input); + if (!input.HasValue()) throw std::runtime_error("actual input"); + SerialParallelFor parallel; + MklPardisoSolver linear; + MklPardisoTangentSolver tangent; + LinearStaticProcedureExecutor oldlinear(parallel, linear); + NonlinearStaticProcedureExecutor oldnonlinear(parallel, tangent); + LinearFrequencyProcedureExecutor oldfrequency; + SolidStaticProcedureExecutor stat(parallel, linear, tangent); + SolidFrequencyProcedureExecutor frequency; + SolidHhtDynamicProcedureExecutor dynamic(parallel, linear, tangent); + NewtonIterationLogWriter log; + if (!log.Open(files.log).IsOk()) throw std::runtime_error("log"); + auto result = ExecutePreparedAnalysis( + input.Value().domain, oldlinear, oldnonlinear, oldfrequency, log, nullptr, + nullptr, &stat, &frequency, &dynamic); + if (!log.Close().IsOk() || !result.HasValue()) + throw std::runtime_error("actual execute"); + return {std::move(input.Value()), std::move(result.Value())}; +} +SolidStateView View(const Actual& actual) { + if (const auto* value = + std::get_if(&actual.outcome)) + return std::cref(value->state); + if (const auto* value = std::get_if(&actual.outcome)) + return std::cref(*value); + return std::cref(std::get(actual.outcome).state); +} +std::vector Bytes(const std::filesystem::path& path) { + std::ifstream stream(path, std::ios::binary); + return {std::istreambuf_iterator(stream), {}}; +} +void Existing(const std::filesystem::path& path) { + std::ofstream out(path, std::ios::binary); + out << "previous final preserved"; +} +std::string Attr(hid_t object, const char* name) { + Hdf5Handle a{H5Aopen(object, name, H5P_DEFAULT), H5Aclose}; + Hdf5Handle type{H5Aget_type(a.Get()), H5Tclose}; + EXPECT_EQ(H5Tget_class(type.Get()), H5T_STRING); + EXPECT_GT(H5Tis_variable_str(type.Get()), 0); + EXPECT_EQ(H5Tget_cset(type.Get()), H5T_CSET_UTF8); + char* pointer = nullptr; + EXPECT_GE(H5Aread(a.Get(), type.Get(), &pointer), 0); + std::string value = pointer ? pointer : ""; + H5free_memory(pointer); + return value; +} +void Compound(hid_t file, const char* path, + const std::vector& names, + const std::vector& kinds, + const std::vector& widths = {}) { + Hdf5Handle dataset{H5Dopen2(file, path, H5P_DEFAULT), H5Dclose}; + ASSERT_GE(dataset.Get(), 0); + Hdf5Handle type{H5Dget_type(dataset.Get()), H5Tclose}; + ASSERT_EQ(H5Tget_class(type.Get()), H5T_COMPOUND); + ASSERT_EQ(H5Tget_nmembers(type.Get()), static_cast(names.size())); + for (std::size_t i = 0; i < names.size(); ++i) { + char* name = H5Tget_member_name(type.Get(), static_cast(i)); + ASSERT_NE(name, nullptr); + EXPECT_EQ(name, names[i]); + H5free_memory(name); + Hdf5Handle field{H5Tget_member_type(type.Get(), static_cast(i)), + H5Tclose}; + if (!widths.empty() && widths[i] > 1) { + ASSERT_EQ(H5Tget_class(field.Get()), H5T_ARRAY); + hsize_t width = 0; + ASSERT_EQ(H5Tget_array_dims2(field.Get(), &width), 1); + EXPECT_EQ(width, widths[i]); + field = Hdf5Handle{H5Tget_super(field.Get()), H5Tclose}; + } + if (kinds[i] == 's') { + EXPECT_EQ(H5Tget_class(field.Get()), H5T_STRING); + EXPECT_GT(H5Tis_variable_str(field.Get()), 0); + EXPECT_EQ(H5Tget_cset(field.Get()), H5T_CSET_UTF8); + } else + EXPECT_GT(H5Tequal(field.Get(), kinds[i] == 'u' ? H5T_STD_U64LE + : kinds[i] == 'b' ? H5T_STD_U8LE + : H5T_IEEE_F64LE), + 0); + } +} + +template +T NumberAttr(hid_t object, const char* name, hid_t file_type, + hid_t native_type) { + Hdf5Handle attribute{H5Aopen(object, name, H5P_DEFAULT), H5Aclose}; + Hdf5Handle type{H5Aget_type(attribute.Get()), H5Tclose}; + Hdf5Handle space{H5Aget_space(attribute.Get()), H5Sclose}; + EXPECT_GT(H5Tequal(type.Get(), file_type), 0); + EXPECT_EQ(H5Sget_simple_extent_type(space.Get()), H5S_SCALAR); + T value{}; + EXPECT_GE(H5Aread(attribute.Get(), native_type, &value), 0); + return value; +} +void DatasetContract(hid_t file, const std::string& path, + const std::vector& expected, + const char* components, const char* units, + const char* location) { + Hdf5Handle dataset{H5Dopen2(file, path.c_str(), H5P_DEFAULT), H5Dclose}; + ASSERT_GE(dataset.Get(), 0); + Hdf5Handle type{H5Dget_type(dataset.Get()), H5Tclose}; + EXPECT_GT(H5Tequal(type.Get(), H5T_IEEE_F64LE), 0); + Hdf5Handle space{H5Dget_space(dataset.Get()), H5Sclose}; + ASSERT_EQ(H5Sget_simple_extent_ndims(space.Get()), + static_cast(expected.size())); + std::vector shape(expected.size()); + H5Sget_simple_extent_dims(space.Get(), shape.data(), nullptr); + EXPECT_EQ(shape, expected); + EXPECT_EQ(Attr(dataset.Get(), "component_names"), components); + EXPECT_EQ(Attr(dataset.Get(), "unit_dimensions"), units); + EXPECT_EQ(Attr(dataset.Get(), "location"), location); + EXPECT_EQ(Attr(dataset.Get(), "coordinate_system"), "global-cartesian"); + const auto marker = path.find("/frames/"); + if (marker != std::string::npos) { + const auto index = + static_cast(std::stoull(path.substr(marker + 8))); + Hdf5Handle metadata{H5Gopen2(file, "/metadata", H5P_DEFAULT), H5Gclose}; + const auto procedure = Attr(metadata.Get(), "analysis_procedure"); + const bool dynamic = procedure.find("hht-") == 0; + const bool nonlinear = procedure == "nonlinear-static-full-newton" || + procedure == "hht-nonlinear-dynamics"; + EXPECT_EQ(Attr(dataset.Get(), "step_name"), "Step-1"); + EXPECT_EQ(NumberAttr(dataset.Get(), "frame_index", + H5T_STD_U64LE, H5T_NATIVE_UINT64), + index); + EXPECT_EQ(NumberAttr(dataset.Get(), "increment_index", + H5T_STD_U64LE, H5T_NATIVE_UINT64), + dynamic ? index + : procedure == "linear-static" ? 1 + : index + 1); + EXPECT_EQ( + NumberAttr(dataset.Get(), dynamic ? "time" : "load_factor", + H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE), + dynamic ? static_cast(index) * .1 + : procedure == "linear-static" ? 1 + : static_cast(index + 1) * .5); + if (path.find("/solid/stress") != std::string::npos) + EXPECT_EQ(Attr(dataset.Get(), "stress_measure"), "cauchy"); + if (path.find("/solid/strain") != std::string::npos) { + EXPECT_EQ(Attr(dataset.Get(), "strain_measure"), + nonlinear ? "green-lagrange" : "infinitesimal"); + EXPECT_EQ(Attr(dataset.Get(), "shear_convention"), "tensor"); + } + if (path.find("/nodal/acceleration") != std::string::npos) + EXPECT_EQ(Attr(dataset.Get(), "acceleration_kind"), "physical"); + if (path.find("/nodal/algorithmic_acceleration") != std::string::npos) + EXPECT_EQ(Attr(dataset.Get(), "acceleration_kind"), "algorithmic"); + if (path.find("/global/strain_energy") != std::string::npos) + EXPECT_EQ(Attr(dataset.Get(), "definition"), + nonlinear ? "sum(W*detJ0)" : "0.5*u^T*K0*u"); + if (path.find("/global/kinetic_energy") != std::string::npos) + EXPECT_EQ(Attr(dataset.Get(), "definition"), "0.5*v^T*M*v"); + if (path.find("/global/external_work") != std::string::npos) + EXPECT_EQ(Attr(dataset.Get(), "definition"), + "sum(0.5*(F_old+F_new)*delta_u)"); + if (path.find("/global/energy_balance_error") != std::string::npos) + EXPECT_EQ(Attr(dataset.Get(), "definition"), "T+U-T0-U0-Wext"); + if (path.find("/global/angular_momentum") != std::string::npos) + EXPECT_EQ(Attr(dataset.Get(), "origin"), "global0"); + } +} +TEST(SolidHdf5ResultsWriter, + AllFiveActualApiOutcomesHaveIndependentOfficialSchema) { + for (unsigned p = 0; p < 5; ++p) { + SCOPED_TRACE(p); + solid_output_test::Files files(p); + auto actual = Execute(files); + ASSERT_TRUE(WriteSolidFile(files.output, actual.input.domain, View(actual), + actual.input.diagnostics) + .IsOk()); + Hdf5Handle file{ + H5Fopen(files.output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), + H5Fclose}; + ASSERT_GE(file.Get(), 0); + Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), + H5Gclose}; + EXPECT_EQ(Attr(metadata.Get(), "analysis_procedure"), + solid_output_test::Procedure(p)); + EXPECT_EQ(Attr(metadata.Get(), "feature_id"), "hex8-structural-analysis"); + EXPECT_EQ(Attr(metadata.Get(), "source_step_name"), "RawSolid"); + EXPECT_EQ(Attr(metadata.Get(), "element_formulation"), "FESA-Hex8"); + EXPECT_EQ(Attr(metadata.Get(), "integration_rule"), "gauss-2x2x2"); + EXPECT_EQ(Attr(metadata.Get(), "mass_formulation"), "consistent-reference"); + EXPECT_EQ(Attr(metadata.Get(), "component_order"), "node-major; UX,UY,UZ"); + Compound( + file.Get(), "/model/nodes", + {"internal_node_id", "instance_name", "source_label", "coordinates"}, + {'u', 's', 's', 'd'}, {1, 1, 1, 3}); + Compound( + file.Get(), "/model/elements", + {"internal_element_id", "instance_name", "source_label", + "source_element_type", "internal_formulation", "node_internal_ids", + "solid_section_internal_id", "material_internal_id"}, + {'u', 's', 's', 's', 's', 'u', 'u', 'u'}, {1, 1, 1, 1, 1, 8, 1, 1}); + Compound(file.Get(), "/model/solid/materials", + {"internal_material_id", "name", "youngs_modulus", "poisson_ratio", + "density_present", "density"}, + {'u', 's', 'd', 'd', 'b', 'd'}); + Compound(file.Get(), "/model/solid/sections", + {"internal_section_id", "source_file", "source_line", + "source_elset", "material_internal_id"}, + {'u', 's', 'u', 's', 'u'}); + Compound(file.Get(), "/model/solid_integration_points", + {"internal_element_id", "point_id", "natural_coordinates", + "reference_coordinates"}, + {'u', 'b', 'd', 'd'}, {1, 1, 3, 3}); + Compound(file.Get(), "/diagnostics", + {"severity", "code", "file", "line", "keyword", "entity_identity", + "message"}, + {'s', 's', 's', 'u', 's', 's', 's'}); + Hdf5Handle mask{ + H5Dopen2(file.Get(), "/model/nodal_constraint_mask", H5P_DEFAULT), + H5Dclose}; + Hdf5Handle masktype{H5Dget_type(mask.Get()), H5Tclose}; + EXPECT_GT(H5Tequal(masktype.Get(), H5T_STD_U8LE), 0); + EXPECT_EQ(Attr(mask.Get(), "component_names"), "UX,UY,UZ"); + if (p == 2) { + EXPECT_EQ(H5Lexists(file.Get(), "/steps/Step-1/frames", H5P_DEFAULT), 0); + Compound(file.Get(), "/steps/Step-1/eigen/modes", + {"mode_id", "eigenvalue", "classification", "cluster_id"}, + {'u', 'd', 'b', 'u'}); + Compound(file.Get(), "/steps/Step-1/eigen/frequencies", + {"mode_id", "frequency"}, {'u', 'd'}); + Compound(file.Get(), "/steps/Step-1/eigen/verification/pairs", + {"sorted_rank", "solver_column", "eigenvalue", "modal_mass", + "rayleigh_quotient", "eta_r", "eta_norm", "eta_energy"}, + {'u', 'u', 'd', 'd', 'd', 'd', 'd', 'd'}); + DatasetContract(file.Get(), "/steps/Step-1/eigen/vectors", {1, 8, 3}, + "UX,UY,UZ", "mass-normalized", "nodal"); + for (const auto& quantity : + std::vector>{ + {"modes", "eigenvalue", "1/time^2"}, + {"frequencies", "frequency", "1/time"}}) { + const auto path = + std::string("/steps/Step-1/eigen/") + std::get<0>(quantity); + Hdf5Handle dataset{H5Dopen2(file.Get(), path.c_str(), H5P_DEFAULT), + H5Dclose}; + EXPECT_EQ(Attr(dataset.Get(), "component_names"), + std::get<1>(quantity)); + EXPECT_EQ(Attr(dataset.Get(), "unit_dimensions"), + std::get<2>(quantity)); + EXPECT_EQ(Attr(dataset.Get(), "location"), "global"); + EXPECT_EQ(Attr(dataset.Get(), "coordinate_system"), "global-cartesian"); + EXPECT_EQ(H5Aexists(dataset.Get(), "frame_index"), 0); + EXPECT_EQ(H5Aexists(dataset.Get(), "step_name"), 0); + } + const auto& state = std::get(actual.outcome); + ASSERT_EQ(state.eigenpair_checks[0].solver_column, 0u); + } else { + const unsigned frames = p == 0 ? 1 : p == 1 ? 2 : 3; + if (p == 1) + Compound(file.Get(), "/steps/Step-1/accepted_frames", + {"frame_index", "increment_index", "load_factor", + "iteration_count", "is_factor_one"}, + {'u', 'u', 'd', 'u', 'b'}); + if (p >= 3) + Compound(file.Get(), "/steps/Step-1/accepted_frames", + {"frame_index", "increment_index", "time", "iteration_count", + "is_initial"}, + {'u', 'u', 'd', 'u', 'b'}); + if (p == 0) + EXPECT_EQ( + H5Lexists(file.Get(), "/steps/Step-1/accepted_frames", H5P_DEFAULT), + 0); + for (unsigned i = 0; i < frames; ++i) { + const auto root = "/steps/Step-1/frames/" + std::to_string(i); + DatasetContract(file.Get(), root + "/nodal/displacement", {8, 3}, + "U1,U2,U3", "length", "nodal"); + DatasetContract(file.Get(), root + "/nodal/reaction", {8, 3}, + "RF1,RF2,RF3", "force", "nodal"); + DatasetContract(file.Get(), root + "/nodal/residual", {8, 3}, + "R1,R2,R3", "force", "nodal"); + DatasetContract(file.Get(), root + "/solid/stress", {1, 8, 6}, + "S11,S22,S33,S12,S13,S23", "force/length^2", + "integration-point"); + DatasetContract(file.Get(), root + "/solid/strain", {1, 8, 6}, + "E11,E22,E33,E12,E13,E23", "dimensionless", + "integration-point"); + DatasetContract(file.Get(), root + "/global/strain_energy", {}, + "scalar", "force*length", "global"); + if (p >= 3) { + DatasetContract(file.Get(), root + "/nodal/velocity", {8, 3}, + "V1,V2,V3", "length/time", "nodal"); + DatasetContract(file.Get(), root + "/nodal/acceleration", {8, 3}, + "A1,A2,A3", "length/time^2", "nodal"); + DatasetContract(file.Get(), root + "/nodal/algorithmic_acceleration", + {8, 3}, "A1,A2,A3", "length/time^2", "nodal"); + for (const auto* energy : + {"kinetic_energy", "external_work", "energy_balance_error"}) + DatasetContract(file.Get(), root + "/global/" + energy, {}, + "scalar", "force*length", "global"); + DatasetContract(file.Get(), root + "/global/linear_momentum", {3}, + "P1,P2,P3", "mass*length/time", "global"); + DatasetContract(file.Get(), root + "/global/angular_momentum", {3}, + "L1,L2,L3", "mass*length^2/time", "global"); + EXPECT_LE( + H5Lexists(file.Get(), (root + "/nodal/displacement_low").c_str(), + H5P_DEFAULT), + 0); + } + } + } + } +} +TEST(SolidHdf5ResultsWriter, + EveryProcedurePreservesPreviousFinalAcrossLifecycleFailures) { + for (unsigned p = 0; p < 5; ++p) { + solid_output_test::Files files(p); + auto actual = Execute(files); + Existing(files.output); + const auto previous = Bytes(files.output); + for (unsigned failure = 0; failure < 5; ++failure) { + SCOPED_TRACE(p); + SCOPED_TRACE(failure); + SolidFileHooks hooks; + if (failure == 0) + hooks.before_write = [](hid_t) { + throw std::runtime_error("injected write"); + }; + if (failure == 1) + hooks.candidate_close = [](hid_t id) -> herr_t { + H5Fclose(id); + return -1; + }; + if (failure == 2) + hooks.before_self_check = [](const auto& path) { + Hdf5Handle file{ + H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), + H5Fclose}; + H5Ldelete(file.Get(), "/model/elements", H5P_DEFAULT); + }; + if (failure == 3) + hooks.reader_close = [](hid_t id) -> herr_t { + H5Fclose(id); + return -1; + }; + if (failure == 4) + hooks.before_finalize = [](const auto& path) { + std::filesystem::remove(path); + }; + const auto status = + WriteSolidFile(files.output, actual.input.domain, View(actual), + actual.input.diagnostics, hooks); + EXPECT_FALSE(status.IsOk()); + EXPECT_EQ(Bytes(files.output), previous); + for (const auto& entry : + std::filesystem::directory_iterator(files.folder)) + EXPECT_EQ(entry.path().filename().string().find("candidate"), + std::string::npos); + } + } +} +TEST(SolidHdf5ResultsWriter, + MalformedWholePairsAndIdentityFailBeforeCandidateCreation) { + solid_output_test::Files files(3); + auto actual = Execute(files); + Existing(files.output); + const auto previous = Bytes(files.output); + auto& state = std::get(actual.outcome).state; + auto& frame = const_cast( + state.AcceptedFrames().back()); + const auto original = frame; + for (unsigned malformed = 0; malformed < 5; ++malformed) { + SCOPED_TRACE(malformed); + frame = original; + if (malformed == 0) frame.low.displacement = Vector{1}; + if (malformed == 1) + frame.low.displacement[21] = std::numeric_limits::quiet_NaN(); + if (malformed == 2) frame.low.displacement[21] = frame.displacement[21]; + if (malformed == 3) frame.nodal_source_ids[0].source_label_text = "wrong"; + if (malformed == 4) frame.low.residual[21] = 1; + EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, View(actual), + actual.input.diagnostics) + .IsOk()); + EXPECT_EQ(Bytes(files.output), previous); + } + frame = original; + const int rounding = std::fegetround(); + ASSERT_EQ(std::fesetround(FE_DOWNWARD), 0); + const auto status = WriteSolidFile(files.output, actual.input.domain, + View(actual), actual.input.diagnostics); + std::fesetround(rounding); + EXPECT_FALSE(status.IsOk()); + EXPECT_EQ(Bytes(files.output), previous); +} +TEST(SolidHdf5ResultsWriter, + ModalColumnsRemainZeroBasedAndRejectInvalidInventory) { + solid_output_test::Files files(2); + auto actual = Execute(files); + auto& state = std::get(actual.outcome); + ASSERT_EQ(state.eigenpair_checks[0].solver_column, 0u); + Existing(files.output); + const auto previous = Bytes(files.output); + const auto original = state; + state.eigenpair_checks[0].solver_column = state.eigenpair_checks.size(); + EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, View(actual), + actual.input.diagnostics) + .IsOk()); + EXPECT_EQ(Bytes(files.output), previous); + state = original; + state.backend.column_ids[0] = state.eigenpair_checks.size(); + state.eigenpair_checks[0].solver_column = state.backend.column_ids[0]; + EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, View(actual), + actual.input.diagnostics) + .IsOk()); +} + +TEST(SolidHdf5ResultsWriter, + ReadbackRejectsMalformedTypesShapesValuesAndSourceIds) { + solid_output_test::Files files(3); + auto actual = Execute(files); + Existing(files.output); + const auto previous = Bytes(files.output); + for (unsigned malformed = 0; malformed < 9; ++malformed) { + SCOPED_TRACE(malformed); + SolidFileHooks hooks; + hooks.before_self_check = [&](const auto& path) { + Hdf5Handle file{H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), + H5Fclose}; + const char* result = "/steps/Step-1/frames/2/nodal/displacement"; + if (malformed == 0) H5Ldelete(file.Get(), result, H5P_DEFAULT); + if (malformed == 1) { + Hdf5Handle dataset{H5Dopen2(file.Get(), result, H5P_DEFAULT), H5Dclose}; + std::vector values(24, 0); + values[0] = std::numeric_limits::quiet_NaN(); + H5Dwrite(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()); + } + if (malformed == 2 || malformed == 3) { + H5Ldelete(file.Get(), result, H5P_DEFAULT); + hsize_t shape[]{malformed == 2 ? 1u : 8u, 3}; + Hdf5Handle space{H5Screate_simple(2, shape, nullptr), H5Sclose}; + Hdf5Handle data{ + H5Dcreate2(file.Get(), result, + malformed == 2 ? H5T_IEEE_F64LE : H5T_IEEE_F32LE, + space.Get(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT), + H5Dclose}; + } + if (malformed == 4) { + Hdf5Handle dataset{ + H5Dopen2(file.Get(), "/model/solid_integration_points", + H5P_DEFAULT), + H5Dclose}; + Hdf5Handle memory{H5Tcreate(H5T_COMPOUND, 1), H5Tclose}; + H5Tinsert(memory.Get(), "point_id", 0, H5T_NATIVE_UINT8); + std::vector values(8, 99); + H5Dwrite(dataset.Get(), memory.Get(), H5S_ALL, H5S_ALL, H5P_DEFAULT, + values.data()); + } + if (malformed == 5) { + Hdf5Handle dataset{H5Dopen2(file.Get(), "/model/elements", H5P_DEFAULT), + H5Dclose}; + Hdf5Handle memory{H5Tcreate(H5T_COMPOUND, sizeof(std::uint64_t)), + H5Tclose}; + H5Tinsert(memory.Get(), "internal_element_id", 0, H5T_NATIVE_UINT64); + const std::uint64_t value = 99; + H5Dwrite(dataset.Get(), memory.Get(), H5S_ALL, H5S_ALL, H5P_DEFAULT, + &value); + } + if (malformed == 6) { + Hdf5Handle dataset{H5Dopen2(file.Get(), result, H5P_DEFAULT), H5Dclose}; + H5Adelete(dataset.Get(), "unit_dimensions"); + WriteStringAttribute(dataset.Get(), "unit_dimensions", "wrong"); + } + if (malformed == 7) { + Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), + H5Gclose}; + H5Adelete(metadata.Get(), "analysis_procedure"); + Hdf5Handle type{H5Tcopy(H5T_C_S1), H5Tclose}; + H5Tset_size(type.Get(), 16); + Hdf5Handle space{H5Screate(H5S_SCALAR), H5Sclose}; + Hdf5Handle attr{ + H5Acreate2(metadata.Get(), "analysis_procedure", type.Get(), + space.Get(), H5P_DEFAULT, H5P_DEFAULT), + H5Aclose}; + } + if (malformed == 8) { + H5Ldelete(file.Get(), "/model/elements", H5P_DEFAULT); + Hdf5Handle type{H5Tcreate(H5T_COMPOUND, 16), H5Tclose}; + H5Tinsert(type.Get(), "material_internal_id", 0, H5T_STD_U64LE); + H5Tinsert(type.Get(), "internal_element_id", 8, H5T_STD_U64LE); + hsize_t one = 1; + Hdf5Handle space{H5Screate_simple(1, &one, nullptr), H5Sclose}; + Hdf5Handle dataset{ + H5Dcreate2(file.Get(), "/model/elements", type.Get(), space.Get(), + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT), + H5Dclose}; + } + }; + EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, View(actual), + actual.input.diagnostics, hooks) + .IsOk()); + EXPECT_EQ(Bytes(files.output), previous); + } +} +TEST(SolidHdf5ResultsWriter, + ReadbackMatchesActualAcceptedValuesAndFullPairObservables) { + for (unsigned p = 0; p < 5; ++p) { + SCOPED_TRACE(p); + solid_output_test::Files files(p); + auto actual = Execute(files); + ASSERT_TRUE(WriteSolidFile(files.output, actual.input.domain, View(actual), + actual.input.diagnostics) + .IsOk()); + Hdf5Handle file{ + H5Fopen(files.output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), + H5Fclose}; + const auto compare = [&](const std::string& path, + const std::vector& expected) { + Hdf5Handle dataset{H5Dopen2(file.Get(), path.c_str(), H5P_DEFAULT), + H5Dclose}; + ASSERT_GE(dataset.Get(), 0); + std::vector values(expected.size()); + ASSERT_GE(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), + 0); + EXPECT_EQ(values, expected); + }; + const auto vector = [](const Vector& input) { + std::vector values; + for (std::size_t i = 0; i < input.Size(); ++i) values.push_back(input[i]); + return values; + }; + Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), + H5Gclose}; + EXPECT_EQ(Attr(metadata.Get(), "source_input_identity"), + SourceInputIdentity(actual.input.domain)); + if (p == 2) { + const auto& state = std::get(actual.outcome); + std::vector values; + for (std::size_t i = 0; i < state.vectors.Rows(); ++i) + for (std::size_t j = 0; j < state.vectors.Columns(); ++j) + values.push_back(state.vectors(i, j)); + compare("/steps/Step-1/eigen/vectors", values); + } else { + const auto frame = [&](const auto& accepted) { + const auto root = + "/steps/Step-1/frames/" + std::to_string(accepted.frame_index); + compare(root + "/nodal/displacement", vector(accepted.displacement)); + compare(root + "/nodal/reaction", vector(accepted.reaction)); + compare(root + "/nodal/residual", vector(accepted.residual)); + std::vector stress, strain; + for (const auto& point : accepted.solid_results) { + stress.insert(stress.end(), point.cauchy_stress.begin(), + point.cauchy_stress.end()); + strain.insert(strain.end(), point.strain.begin(), point.strain.end()); + } + compare(root + "/solid/stress", stress); + compare(root + "/solid/strain", strain); + }; + if (p < 2) + for (const auto& accepted : + std::get(actual.outcome) + .state.AcceptedFrames()) + frame(accepted); + else { + bool nonzero_low = false; + for (const auto& accepted : + std::get(actual.outcome) + .state.AcceptedFrames()) { + frame(accepted); + const auto root = + "/steps/Step-1/frames/" + std::to_string(accepted.frame_index); + compare(root + "/nodal/velocity", vector(accepted.velocity)); + compare(root + "/nodal/acceleration", vector(accepted.acceleration)); + compare(root + "/nodal/algorithmic_acceleration", + vector(accepted.algorithmic_acceleration)); + compare(root + "/global/strain_energy", {accepted.strain_energy}); + compare(root + "/global/kinetic_energy", {accepted.kinetic_energy}); + compare(root + "/global/external_work", {accepted.external_work}); + compare(root + "/global/energy_balance_error", + {accepted.energy_balance_error}); + for (std::size_t i = 0; i < accepted.low.displacement.Size(); ++i) + nonzero_low = nonzero_low || accepted.low.displacement[i] != 0; + } + EXPECT_TRUE(nonzero_low); + } + } + } +} +class FailedTangent final : public TangentSolver { + public: + Status FactorizeFresh(const SparseMatrix&) override { + return Status::Failure(FailureCategory::kSolver, + {{Severity::kError, + "injected-tangent-failure", + {}, + "", + {}, + "test backend failed"}}); + } + Status SolveCorrection(const Vector&, Vector&) const override { + return Status::Failure(FailureCategory::kSolver, {}); + } +}; +TEST(SolidPreparedAnalysis, + FailureRecordsAreLoggedAndLogFailureTakesPrecedence) { + for (unsigned p : {1u, 3u, 4u}) { + solid_output_test::Files files(p); + auto input = PrepareAnalysisInput(files.input); + ASSERT_TRUE(input.HasValue()); + SerialParallelFor parallel; + MklPardisoSolver linear; + FailedTangent tangent; + LinearStaticProcedureExecutor oldlinear(parallel, linear); + NonlinearStaticProcedureExecutor oldnonlinear(parallel, tangent); + LinearFrequencyProcedureExecutor oldfrequency; + SolidStaticProcedureExecutor stat(parallel, linear, tangent); + SolidFrequencyProcedureExecutor frequency; + SolidHhtDynamicProcedureExecutor dynamic(parallel, linear, tangent); + for (bool closed : {false, true}) { + SCOPED_TRACE(p); + SCOPED_TRACE(closed); + NewtonIterationLogWriter log; + ASSERT_TRUE(log.Open(files.log).IsOk()); + if (closed) ASSERT_TRUE(log.Close().IsOk()); + const auto result = ExecutePreparedAnalysis( + input.Value().domain, oldlinear, oldnonlinear, oldfrequency, log, + nullptr, nullptr, &stat, &frequency, &dynamic); + ASSERT_FALSE(result.HasValue()); + EXPECT_EQ(result.GetStatus().Category(), + closed ? FailureCategory::kOutput : FailureCategory::kSolver); + if (!closed) { + ASSERT_TRUE(log.Close().IsOk()); + const auto bytes = Bytes(files.log); + const std::string text(bytes.begin(), bytes.end()); + EXPECT_NE(text.find("injected-tangent-failure"), std::string::npos); + EXPECT_EQ(text.find("rotation"), std::string::npos); + EXPECT_NE(text.find("stage="), std::string::npos); + } + } + } +} +TEST(SolidPreparedAnalysis, MissingSolidAdapterFailsWithoutLegacyDispatch) { + solid_output_test::Files files(0); + auto input = PrepareAnalysisInput(files.input); + ASSERT_TRUE(input.HasValue()); + SerialParallelFor parallel; + MklPardisoSolver linear; + MklPardisoTangentSolver tangent; + LinearStaticProcedureExecutor oldlinear(parallel, linear); + NonlinearStaticProcedureExecutor oldnonlinear(parallel, tangent); + LinearFrequencyProcedureExecutor oldfrequency; + NewtonIterationLogWriter log; + ASSERT_TRUE(log.Open(files.log).IsOk()); + const auto result = ExecutePreparedAnalysis(input.Value().domain, oldlinear, + oldnonlinear, oldfrequency, log); + ASSERT_FALSE(result.HasValue()); + ASSERT_FALSE(result.GetStatus().Diagnostics().empty()); + EXPECT_EQ(result.GetStatus().Diagnostics()[0].code, "missing-solid-executor"); + EXPECT_TRUE(log.Close().IsOk()); +} + +TEST(SolidHdf5ResultsWriter, UnsupportedEnvironmentAfterCallbackCannotPublish) { + solid_output_test::Files files(3); + auto actual = Execute(files); + Existing(files.output); + const auto previous = Bytes(files.output); + for (unsigned callback = 0; callback < 2; ++callback) { + for (unsigned control : {0x8000u, 0x40u}) { + SCOPED_TRACE(callback); + SCOPED_TRACE(control); + struct Restore { + unsigned original = _mm_getcsr(); + ~Restore() { _mm_setcsr(original); } + } restore; + SolidFileHooks hooks; + if (callback == 0) + hooks.before_write = [&](hid_t) { + _mm_setcsr(restore.original | control); + }; + else + hooks.before_finalize = [&](const auto&) { + _mm_setcsr(restore.original | control); + }; + EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, + View(actual), actual.input.diagnostics, hooks) + .IsOk()); + EXPECT_EQ(Bytes(files.output), previous); + } + } +} +TEST(SolidPreparedAnalysis, MissingWriterPreservesFinalForEveryTypedOutcome) { + for (unsigned p = 0; p < 5; ++p) { + solid_output_test::Files files(p); + auto actual = Execute(files); + Existing(files.output); + const auto previous = Bytes(files.output); + const auto paths = + DeriveSingleInputOutputPaths({files.input.string()}, files.folder); + ASSERT_TRUE(paths.HasValue()); + Hdf5ResultsWriter linear; + NonlinearHdf5ResultsWriter nonlinear; + EigenHdf5ResultsWriter eigen; + const auto status = FinalizePreparedAnalysisResult( + paths.Value(), actual.input.domain, actual.outcome, + actual.input.diagnostics, linear, nonlinear, eigen); + ASSERT_FALSE(status.IsOk()); + ASSERT_FALSE(status.Diagnostics().empty()); + EXPECT_EQ(status.Diagnostics()[0].code, "missing-solid-writer"); + EXPECT_EQ(Bytes(files.output), previous); + } +} +TEST(SolidNewtonIterationLog, + InvalidScalarEvidencePermanentlyFailsLogLifecycle) { + for (unsigned failure = 0; failure < 3; ++failure) { + solid_output_test::Files files(3); + NewtonIterationLogWriter log; + ASSERT_TRUE(log.Open(files.log).IsOk()); + SolidDynamicIterationRecord record; + record.stage = "failed"; + record.outcome = "failed"; + if (failure == 0) + record.evidence.correction_metrics[2] = + std::numeric_limits::infinity(); + if (failure == 1) record.evidence.residual_scale = -1; + if (failure == 2) + record.evidence.worst_weighted_terms[0] = + std::numeric_limits::quiet_NaN(); + EXPECT_FALSE(log.Write(record).IsOk()); + EXPECT_FALSE(log.Close().IsOk()); + } +} +TEST(SolidHdf5ResultsWriter, AbsentStaticDensityIsAnExplicitPlaceholder) { + solid_output_test::Files files(0); + auto deck = solid_output_test::Deck(0); + const auto position = deck.find("*Density\n2\n"); + ASSERT_NE(position, std::string::npos); + deck.erase(position, std::string("*Density\n2\n").size()); + { + std::ofstream out(files.input); + out << deck; + } + auto actual = Execute(files); + ASSERT_TRUE(WriteSolidFile(files.output, actual.input.domain, View(actual), + actual.input.diagnostics) + .IsOk()); + Hdf5Handle file{ + H5Fopen(files.output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), + H5Fclose}; + Hdf5Handle dataset{ + H5Dopen2(file.Get(), "/model/solid/materials", H5P_DEFAULT), H5Dclose}; + struct Density { + std::uint8_t present; + double rho; + } row{}; + Hdf5Handle type{H5Tcreate(H5T_COMPOUND, sizeof(row)), H5Tclose}; + H5Tinsert(type.Get(), "density_present", HOFFSET(Density, present), + H5T_NATIVE_UINT8); + H5Tinsert(type.Get(), "density", HOFFSET(Density, rho), H5T_NATIVE_DOUBLE); + ASSERT_GE( + H5Dread(dataset.Get(), type.Get(), H5S_ALL, H5S_ALL, H5P_DEFAULT, &row), + 0); + EXPECT_EQ(row.present, 0); + EXPECT_EQ(row.rho, 0); +} + +TEST(SolidHdf5ResultsWriter, + ModalPhysicalTableAttributeCorruptionPreservesFinal) { + solid_output_test::Files files(2); + auto actual = Execute(files); + Existing(files.output); + const auto previous = Bytes(files.output); + for (const auto* quantity : {"modes", "frequencies"}) { + SolidFileHooks hooks; + hooks.before_self_check = [&](const auto& path) { + Hdf5Handle file{H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), + H5Fclose}; + const auto dataset_path = std::string("/steps/Step-1/eigen/") + quantity; + Hdf5Handle dataset{ + H5Dopen2(file.Get(), dataset_path.c_str(), H5P_DEFAULT), H5Dclose}; + H5Adelete(dataset.Get(), "unit_dimensions"); + WriteStringAttribute(dataset.Get(), "unit_dimensions", "wrong"); + }; + EXPECT_FALSE(WriteSolidFile(files.output, actual.input.domain, View(actual), + actual.input.diagnostics, hooks) + .IsOk()); + EXPECT_EQ(Bytes(files.output), previous); + } +} + +TEST(SolidHdf5ResultsWriter, + UnselectedBackendPairsCannotPublishInvalidSpectrumOrDimension) { + solid_output_test::Files files(2); + auto deck = solid_output_test::Deck(2); + const auto position = deck.find("Tip,2,3\n"); + ASSERT_NE(position, std::string::npos); + deck.erase(position, std::string("Tip,2,3\n").size()); + { + std::ofstream out(files.input); + out << deck; + } + auto actual = Execute(files); + auto& state = std::get(actual.outcome); + ASSERT_EQ(state.modes.size(), 1u); + ASSERT_EQ(state.free_dof_count, 3u); + ASSERT_GT(state.eigenpair_checks.size(), state.modes.size()); + ASSERT_LE(state.eigenpair_checks.size(), state.free_dof_count); + const auto original = state; + Existing(files.output); + const auto previous = Bytes(files.output); + for (unsigned failure = 0; failure < 3; ++failure) { + SCOPED_TRACE(failure); + state = original; + if (failure == 0) + state.eigenpair_checks[1].eigenvalue = -1e-6 * state.classification_scale; + if (failure == 1) + state.eigenpair_checks[1].eigenvalue = + state.eigenpair_checks[0].eigenvalue * .5; + if (failure < 2) + state.backend.eigenvalues[1] = state.eigenpair_checks[1].eigenvalue; + if (failure == 2) { + while (state.eigenpair_checks.size() <= state.free_dof_count) { + auto pair = state.eigenpair_checks.back(); + pair.sorted_rank = state.eigenpair_checks.size() + 1; + pair.solver_column = state.eigenpair_checks.size(); + state.eigenpair_checks.push_back(pair); + state.backend.eigenvalues.push_back(pair.eigenvalue); + state.backend.column_ids.push_back(pair.solver_column); + } + } + const auto status = WriteSolidFile(files.output, actual.input.domain, + View(actual), actual.input.diagnostics); + ASSERT_FALSE(status.IsOk()); + ASSERT_FALSE(status.Diagnostics().empty()); + const auto& message = status.Diagnostics()[0].message; + EXPECT_NE(message.find(failure == 0 ? "Significant negative" + : failure == 1 ? "not sorted" + : "inventory"), + std::string::npos); + EXPECT_EQ(Bytes(files.output), previous); + } +} +} // namespace diff --git a/tests/unit/io/hdf5/solid_output_test_fixture.h b/tests/unit/io/hdf5/solid_output_test_fixture.h new file mode 100644 index 0000000..b5dbf13 --- /dev/null +++ b/tests/unit/io/hdf5/solid_output_test_fixture.h @@ -0,0 +1,58 @@ +#ifndef FESA_TESTS_SOLID_OUTPUT_FIXTURE_H_ +#define FESA_TESTS_SOLID_OUTPUT_FIXTURE_H_ +#include +#include +#include +#include +#include +namespace solid_output_test { +inline std::string Deck(unsigned procedure) { + std::string s = + "*Part,name=P\n*Node\n10,0,0,0\n20,1,0,0\n30,1,1,0\n40,0,1,0\n" + "50,0,0,1\n60,1,0,1\n70,1,1,1\n80,0,1,1\n" + "*Element,type=C3D8,elset=Solid\n99,10,20,30,40,50,60,70,80\n" + "*Solid Section,elset=Solid,material=M\n,\n*End Part\n" + "*Assembly,name=A\n*Instance,name=I,part=P\n*End Instance\n" + "*Nset,nset=Fixed,instance=I\n10,20,30,40,50,60,70\n" + "*Nset,nset=Tip,instance=I\n80\n*End Assembly\n" + "*Material,name=M\n*Elastic\n100,.25\n*Density\n2\n" + "*Boundary\nFixed,1,6\nTip,2,3\n"; + s += "*Step,name=RawSolid,inc=2,nlgeom="; + s += (procedure == 1 || procedure == 4) ? "YES\n" : "NO\n"; + if (procedure == 0) s += "*Static\n1,1,1,1\n"; + if (procedure == 1) s += "*Static,Direct\n.5,1\n"; + if (procedure == 2) + s += "*Frequency,eigensolver=Lanczos,normalization=mass\n1\n"; + if (procedure >= 3) s += "*Dynamic,Direct\n.1,.2\n"; + if (procedure != 2) s += "*Cload\nTip,1,.001\n"; + return s + "*End Step\n"; +} +struct Files { + std::filesystem::path folder, input, output, log; + explicit Files(unsigned procedure) { + static std::atomic sequence{0}; + folder = std::filesystem::temp_directory_path() / + ("fesa-solid-output-" + + std::to_string( + std::chrono::steady_clock::now().time_since_epoch().count()) + + "-" + std::to_string(sequence.fetch_add(1))); + std::filesystem::create_directory(folder); + input = folder / "Solid.inp"; + output = folder / "Solid_results.h5"; + log = folder / "Solid_log.txt"; + std::ofstream out(input); + out << Deck(procedure); + } + ~Files() { + std::error_code ignored; + std::filesystem::remove_all(folder, ignored); + } +}; +inline const char* Procedure(unsigned p) { + static const char* names[] = {"linear-static", "nonlinear-static-full-newton", + "linear-frequency", "hht-linear-dynamics", + "hht-nonlinear-dynamics"}; + return names[p]; +} +} // namespace solid_output_test +#endif