diff --git a/src/fesa/io/hdf5/hdf5_results_writer.cpp b/src/fesa/io/hdf5/hdf5_results_writer.cpp index b3b58b5..54aea71 100644 --- a/src/fesa/io/hdf5/hdf5_results_writer.cpp +++ b/src/fesa/io/hdf5/hdf5_results_writer.cpp @@ -3,7 +3,9 @@ #include "fesa/io/hdf5/hdf5_results_writer.hpp" +#include "fesa/analysis/analysis_model.hpp" #include "fesa/build_info.hpp" +#include "fesa/fem/dof_manager.hpp" #include @@ -28,6 +30,11 @@ constexpr std::size_t kEndpointCount = 2U; constexpr std::size_t kGaussPointCount = 2U; constexpr std::size_t kEndActionComponentCount = 6U; constexpr std::size_t kGeneralizedComponentCount = 4U; +constexpr std::size_t kShellNodeCount = 4U; +constexpr std::size_t kShellLocationCount = 4U; +constexpr std::size_t kShellGeneralizedComponentCount = 8U; +constexpr std::size_t kShellSectionPositionCount = 3U; +constexpr std::size_t kShellStressComponentCount = 3U; constexpr const char* kStepName = "Step-1"; constexpr std::size_t kFrameIndex = 0U; constexpr const char* kStepRoot = "/steps/Step-1/frames/0"; @@ -217,6 +224,46 @@ bool sameIdentity(const SourceEntityId& left, const SourceEntityId& right) { using AxisSet = std::array; +struct WriterModelData { + std::vector beamLocalAxes; + std::vector constraintMask; + std::vector prescribedDisplacement; +}; + +bool isShellDomain(const Domain& domain) noexcept { + return !domain.shellElements().empty(); +} + +const char* shellSourceTypeName(const ShellSourceElementType type) { + return type == ShellSourceElementType::s4 ? "S4" : "S4R"; +} + +bool isOrthonormalRightHanded( + const std::array, 3>& frame) { + constexpr double tolerance = 1.0e-12; + const auto dot = [](const std::array& left, + const std::array& right) { + return left[0U] * right[0U] + left[1U] * right[1U] + + left[2U] * right[2U]; + }; + for (const auto& axis : frame) { + if (!isFinite(axis) || std::abs(dot(axis, axis) - 1.0) > tolerance) { + return false; + } + } + if (std::abs(dot(frame[0U], frame[1U])) > tolerance || + std::abs(dot(frame[0U], frame[2U])) > tolerance || + std::abs(dot(frame[1U], frame[2U])) > tolerance) { + return false; + } + const std::array cross = { + frame[0U][1U] * frame[1U][2U] - frame[0U][2U] * frame[1U][1U], + frame[0U][2U] * frame[1U][0U] - frame[0U][0U] * frame[1U][2U], + frame[0U][0U] * frame[1U][1U] - frame[0U][1U] * frame[1U][0U]}; + return dot(cross, frame[2U]) > 0.0 && + std::abs(dot(cross, frame[2U]) - 1.0) <= tolerance; +} + bool computeLocalAxes( const Domain& domain, const EulerBeam3DDefinition& element, AxisSet& axes) { if (element.nodeIndices[0U] >= domain.nodes().size() || @@ -270,12 +317,142 @@ bool sizeProductFits( return true; } +Status validateShellWriterInput( + const Domain& domain, const AnalysisState& state) { + if (!domain.elements().empty()) { + return outputFailure( + "invalid-result-identity", + "Schema v0 does not combine B33 and FESA-MITC4 element inventories."); + } + for (const auto& material : domain.materials()) { + if (material.name.empty() || !isValidUtf8(material.name) || + !std::isfinite(material.youngsModulus) || + !std::isfinite(material.poissonRatio) || + !(material.youngsModulus > 0.0)) { + return outputFailure( + "invalid-result-identity", + "Every shell material requires a UTF-8 name and finite constitutive data."); + } + } + for (const auto& section : domain.shellSections()) { + if (section.name.empty() || !isValidUtf8(section.name) || + section.materialIndex >= domain.materials().size() || + !std::isfinite(section.thickness) || !(section.thickness > 0.0)) { + return outputFailure( + "invalid-result-identity", + "Every shell section requires stable material identity and positive finite thickness."); + } + } + for (const auto& element : domain.shellElements()) { + if ((element.sourceType != ShellSourceElementType::s4 && + element.sourceType != ShellSourceElementType::s4r) || + element.sourceId.sourceLabel <= 0 || + element.sourceId.sourceLabelText.empty() || + !isValidUtf8(element.sourceId.instanceName) || + !isValidUtf8(element.sourceId.sourceLabelText) || + element.materialIndex >= domain.materials().size() || + element.sectionIndex >= domain.shellSections().size() || + domain.shellSections()[element.sectionIndex].materialIndex != + element.materialIndex) { + return outputFailure( + "invalid-result-identity", + "Every shell element requires stable source, material, and section identity."); + } + std::array sortedNodes = element.nodeIndices; + std::sort(sortedNodes.begin(), sortedNodes.end()); + if (sortedNodes.back() >= domain.nodes().size() || + std::adjacent_find(sortedNodes.begin(), sortedNodes.end()) != + sortedNodes.end()) { + return outputFailure( + "invalid-result-identity", + "Every shell element requires four distinct valid node identities."); + } + } + if (domain.shellNodeInitialFrames().size() != domain.nodes().size()) { + return outputFailure( + "invalid-result-identity", + "Shell output requires one initial director/frame per source-ordered node."); + } + for (std::size_t node = 0U; + node < domain.shellNodeInitialFrames().size(); + ++node) { + const auto& source = domain.shellNodeInitialFrames()[node]; + const std::array, 3> frame{ + source.tangentA, source.tangentB, source.director}; + if (source.nodeIndex != node || !isOrthonormalRightHanded(frame)) { + return outputFailure( + "invalid-result-identity", + "Shell initial frames must be finite, orthonormal, right-handed, and node ordered."); + } + } + + std::size_t expectedRows = 0U; + if (!sizeProductFits( + domain.shellElements().size(), kShellLocationCount, expectedRows) || + state.shellResults().size() != expectedRows) { + return outputFailure( + "invalid-result-rows", + "Shell output requires exactly GP1 through GP4 for every shell element."); + } + const double gauss = 1.0 / std::sqrt(3.0); + const std::array locations{ + ShellMidsurfaceLocation::gp1, + ShellMidsurfaceLocation::gp2, + ShellMidsurfaceLocation::gp3, + ShellMidsurfaceLocation::gp4}; + const std::array, kShellLocationCount> coordinates{{ + {-gauss, -gauss}, + {gauss, -gauss}, + {gauss, gauss}, + {-gauss, gauss}}}; + const std::array positions{ + ShellSectionPosition::bottom, + ShellSectionPosition::middle, + ShellSectionPosition::top}; + constexpr std::array zeta{-1.0, 0.0, 1.0}; + for (std::size_t rowIndex = 0U; + rowIndex < state.shellResults().size(); + ++rowIndex) { + const auto& row = state.shellResults()[rowIndex]; + const std::size_t element = rowIndex / kShellLocationCount; + const std::size_t location = rowIndex % kShellLocationCount; + if (row.element != element || row.location != locations[location] || + row.naturalCoordinates != coordinates[location] || + !isOrthonormalRightHanded(row.localFrame) || + !isFinite(row.generalizedStrain) || + !isFinite(row.sectionResultant)) { + return outputFailure( + "invalid-result-rows", + "Shell result rows must be finite and preserve element/GP/frame identity."); + } + for (std::size_t position = 0U; + position < kShellSectionPositionCount; + ++position) { + if (row.stress[position].position != positions[position] || + row.stress[position].zeta != zeta[position] || + !isFinite(row.stress[position].components)) { + return outputFailure( + "invalid-result-rows", + "Shell stress rows must preserve BOTTOM, MIDDLE, TOP identity."); + } + } + } + if (!std::isfinite(state.physicalStrainEnergy()) || + !isFinite(state.equilibrium()) || + !isFinite(state.verificationMetrics())) { + return outputFailure( + "invalid-result-rows", + "Shell energy, equilibrium, and verification metrics must be finite."); + } + return Status::ok(); +} + Status validateWriterInput( const std::filesystem::path& outputPath, const Domain& domain, const AnalysisState& state, const std::vector& diagnostics, - std::vector& localAxes) { + WriterModelData& modelData) { if (outputPath.empty() || outputPath.filename().empty()) { return outputFailure( "invalid-output-path", "The HDF5 output path must name a file."); @@ -287,6 +464,18 @@ Status validateWriterInput( "Schema v0 requires literal Step-1 and frame index 0."); } + const bool shell = isShellDomain(domain); + if (shell) { + const Status shellValidation = validateShellWriterInput(domain, state); + if (!shellValidation.isOk()) { + return shellValidation; + } + } else if (!state.shellResults().empty()) { + return outputFailure( + "invalid-result-rows", + "Beam output cannot contain shell recovery rows."); + } + std::size_t fullDofCount = 0U; if (!sizeProductFits(domain.nodes().size(), kDofsPerNode, fullDofCount)) { return outputFailure( @@ -331,8 +520,8 @@ Status validateWriterInput( } } - localAxes.clear(); - localAxes.reserve(domain.elements().size()); + modelData.beamLocalAxes.clear(); + modelData.beamLocalAxes.reserve(domain.elements().size()); for (const EulerBeam3DDefinition& element : domain.elements()) { AxisSet axes{}; if (element.sourceId.sourceLabel <= 0 || @@ -346,7 +535,7 @@ Status validateWriterInput( "invalid-result-identity", "Every element requires valid source, connectivity, property, and local-axis identity."); } - localAxes.push_back(axes); + modelData.beamLocalAxes.push_back(axes); } std::size_t endpointCount = 0U; @@ -444,6 +633,43 @@ Status validateWriterInput( "Diagnostic text must be valid UTF-8."); } } + + + auto analysisModelResult = AnalysisModel::create(domain); + if (!analysisModelResult.hasValue()) { + return outputFailure( + "invalid-result-state", + "The HDF5 writer could not reconstruct the active model view."); + } + const AnalysisModel analysisModel = + std::move(analysisModelResult.value()); + auto dofResult = DofManager::create(analysisModel); + if (!dofResult.hasValue()) { + return outputFailure( + "invalid-result-state", + "The HDF5 writer could not reconstruct stable constraint identity."); + } + const DofManager dofs = std::move(dofResult.value()); + modelData.constraintMask.assign(fullDofCount, 0U); + modelData.prescribedDisplacement.assign(fullDofCount, 0.0); + if (dofs.constrainedDofs().size() != dofs.prescribedValues().size()) { + return outputFailure( + "invalid-result-state", + "Constraint identities and prescribed values have inconsistent sizes."); + } + for (std::size_t index = 0U; + index < dofs.constrainedDofs().size(); + ++index) { + const std::size_t fullDof = dofs.constrainedDofs()[index]; + const double prescribed = dofs.prescribedValues()[index]; + if (fullDof >= fullDofCount || !std::isfinite(prescribed)) { + return outputFailure( + "invalid-result-state", + "Constraint identity and prescribed values must be finite and in range."); + } + modelData.constraintMask[fullDof] = 1U; + modelData.prescribedDisplacement[fullDof] = prescribed; + } return Status::ok(); } @@ -545,16 +771,12 @@ void writeResultAttributes( writeUint64Attribute(dataset, "frame_index", 0U); } -void writeDoubleDataset( +Hdf5Handle writeDoubleValues( const hid_t file, const std::string& path, const std::vector& dimensions, const double* values, - const std::size_t valueCount, - const std::string& componentNames, - const std::string& componentUnits, - const std::string& coordinateSystem, - const std::string& location) { + const std::size_t valueCount) { auto space = createDatasetSpace(dimensions); Hdf5Handle dataset{ requireHdf5Id( @@ -570,10 +792,69 @@ void writeDoubleDataset( H5P_DEFAULT, values), "Unable to write a floating-point HDF5 dataset."); } + return dataset; +} + +void writeDoubleDataset( + const hid_t file, + const std::string& path, + const std::vector& dimensions, + const double* values, + const std::size_t valueCount, + const std::string& componentNames, + const std::string& componentUnits, + const std::string& coordinateSystem, + const std::string& location) { + auto dataset = writeDoubleValues( + file, path, dimensions, values, valueCount); writeResultAttributes( dataset.get(), componentNames, componentUnits, coordinateSystem, location); } +void writeModelDoubleDataset( + const hid_t file, + const std::string& path, + const std::vector& dimensions, + const double* values, + const std::size_t valueCount, + const std::string& componentNames, + const std::string& componentUnits, + const std::string& coordinateSystem, + const std::string& location) { + auto dataset = writeDoubleValues( + file, path, dimensions, values, valueCount); + writeStringAttribute(dataset.get(), "component_names", componentNames); + writeStringAttribute( + dataset.get(), "component_unit_dimensions", componentUnits); + writeStringAttribute(dataset.get(), "coordinate_system", coordinateSystem); + writeStringAttribute(dataset.get(), "location", location); +} + +void writeUint8Dataset( + const hid_t file, + const std::string& path, + const std::vector& dimensions, + const std::uint8_t* values, + const std::size_t valueCount) { + auto space = createDatasetSpace(dimensions); + Hdf5Handle dataset{ + requireHdf5Id( + H5Dcreate2( + file, path.c_str(), H5T_STD_U8LE, space.get(), + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT), + "Unable to create a uint8 HDF5 dataset."), + H5Dclose}; + if (valueCount != 0U) { + requireHdf5( + H5Dwrite( + dataset.get(), H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values), + "Unable to write a uint8 HDF5 dataset."); + } + writeStringAttribute(dataset.get(), "component_names", "UX,UY,UZ,URX,URY,URZ"); + writeStringAttribute(dataset.get(), "value_meaning", "0=free,1=constrained"); +} + struct NodeWriteRow { std::uint64_t internalNodeId; const char* instanceName; @@ -589,6 +870,33 @@ struct ElementWriteRow { double localAxes[9]; }; +struct ShellElementWriteRow { + std::uint64_t internalElementId; + const char* instanceName; + const char* sourceLabel; + const char* sourceElementType; + const char* internalFormulation; + std::uint64_t nodeInternalIds[4]; + std::uint64_t shellSectionInternalId; + std::uint64_t materialInternalId; +}; + +struct ShellMaterialWriteRow { + std::uint64_t internalMaterialId; + const char* name; + double youngsModulus; + double poissonRatio; +}; + +struct ShellSectionWriteRow { + std::uint64_t internalSectionId; + const char* sourceFile; + std::uint64_t sourceLine; + const char* sourceElset; + std::uint64_t materialInternalId; + double thickness; +}; + struct StressWriteRow { std::uint64_t internalElementId; std::uint64_t gaussPointIndex; @@ -639,20 +947,35 @@ Hdf5Handle writeCompoundDataset( void writeMetadata(const hid_t file, const Domain& domain) { auto metadata = createGroup(file, "/metadata"); writeUint64Attribute(metadata.get(), "schema_version", 0U); - writeStringAttribute( - metadata.get(), "feature_id", "linear-static-3d-euler-beam"); writeStringAttribute( metadata.get(), "solver_version", std::string{solverVersion()}); writeStringAttribute( metadata.get(), "source_input_identity", sourceInputIdentity(domain)); writeStringAttribute( metadata.get(), "unit_system_label", "user-consistent-unspecified"); - writeStringAttribute( - metadata.get(), - "coordinate_convention", - "global-cartesian; beam-local=(t,n1,t-cross-n1)"); - writeStringAttribute( - metadata.get(), "element_formulation", "B33-3D-Euler-Bernoulli"); + if (isShellDomain(domain)) { + writeStringAttribute( + metadata.get(), "feature_id", "linear-static-mitc4-shell"); + writeStringAttribute( + metadata.get(), + "coordinate_convention", + "global-cartesian; shell-local=(e1,e2,e3); positive-thickness=+zeta"); + writeStringAttribute( + metadata.get(), "internal_formulation", "FESA-MITC4"); + writeStringAttribute( + metadata.get(), + "integration_rule", + "2x2x2-gauss; mitc4-edge-midpoint-shear"); + } else { + writeStringAttribute( + metadata.get(), "feature_id", "linear-static-3d-euler-beam"); + writeStringAttribute( + metadata.get(), + "coordinate_convention", + "global-cartesian; beam-local=(t,n1,t-cross-n1)"); + writeStringAttribute( + metadata.get(), "element_formulation", "B33-3D-Euler-Bernoulli"); + } writeStringAttribute(metadata.get(), "step_name", kStepName); writeUint64Attribute(metadata.get(), "frame_index", 0U); } @@ -731,7 +1054,7 @@ void writeNodes(const hid_t file, const Domain& domain) { writeStringAttribute(dataset.get(), "units_label", "length"); } -void writeElements( +void writeBeamElements( const hid_t file, const Domain& domain, const std::vector& axes) { std::vector rows; rows.reserve(domain.elements().size()); @@ -816,6 +1139,279 @@ void writeElements( dataset.get(), "formulation", "B33-3D-Euler-Bernoulli"); } +void writeShellElements(const hid_t file, const Domain& domain) { + std::vector rows; + rows.reserve(domain.shellElements().size()); + for (std::size_t index = 0U; index < domain.shellElements().size(); ++index) { + const auto& element = domain.shellElements()[index]; + rows.push_back({ + static_cast(index), + element.sourceId.instanceName.c_str(), + element.sourceId.sourceLabelText.c_str(), + shellSourceTypeName(element.sourceType), + kMitc4InternalFormulation.data(), + {static_cast(element.nodeIndices[0U]), + static_cast(element.nodeIndices[1U]), + static_cast(element.nodeIndices[2U]), + static_cast(element.nodeIndices[3U])}, + static_cast(element.sectionIndex), + static_cast(element.materialIndex)}); + } + + auto stringType = makeUtf8StringType(); + const hsize_t nodeDimensions[] = {kShellNodeCount}; + Hdf5Handle fileNodes{ + requireHdf5Id( + H5Tarray_create2(H5T_STD_U64LE, 1, nodeDimensions), + "Unable to create shell connectivity file type."), + H5Tclose}; + Hdf5Handle memoryNodes{ + requireHdf5Id( + H5Tarray_create2(H5T_NATIVE_UINT64, 1, nodeDimensions), + "Unable to create shell connectivity memory type."), + H5Tclose}; + Hdf5Handle fileType{ + requireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(ShellElementWriteRow)), + "Unable to create shell element file type."), + H5Tclose}; + Hdf5Handle memoryType{ + requireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(ShellElementWriteRow)), + "Unable to create shell element memory type."), + H5Tclose}; + const auto insertFields = [&](const hid_t type, + const hid_t integerType, + const hid_t nodeType) { + requireHdf5( + H5Tinsert(type, "internal_element_id", + HOFFSET(ShellElementWriteRow, internalElementId), integerType), + "Unable to define shell element ID field."); + requireHdf5( + H5Tinsert(type, "instance_name", + HOFFSET(ShellElementWriteRow, instanceName), stringType.get()), + "Unable to define shell element instance field."); + requireHdf5( + H5Tinsert(type, "source_label", + HOFFSET(ShellElementWriteRow, sourceLabel), stringType.get()), + "Unable to define shell element source-label field."); + requireHdf5( + H5Tinsert(type, "source_element_type", + HOFFSET(ShellElementWriteRow, sourceElementType), stringType.get()), + "Unable to define shell source-type field."); + requireHdf5( + H5Tinsert(type, "internal_formulation", + HOFFSET(ShellElementWriteRow, internalFormulation), stringType.get()), + "Unable to define shell formulation field."); + requireHdf5( + H5Tinsert(type, "node_internal_ids", + HOFFSET(ShellElementWriteRow, nodeInternalIds), nodeType), + "Unable to define shell connectivity field."); + requireHdf5( + H5Tinsert(type, "shell_section_internal_id", + HOFFSET(ShellElementWriteRow, shellSectionInternalId), integerType), + "Unable to define shell section ID field."); + requireHdf5( + H5Tinsert(type, "material_internal_id", + HOFFSET(ShellElementWriteRow, materialInternalId), integerType), + "Unable to define shell material ID field."); + }; + insertFields(fileType.get(), H5T_STD_U64LE, fileNodes.get()); + insertFields(memoryType.get(), H5T_NATIVE_UINT64, memoryNodes.get()); + auto dataset = writeCompoundDataset( + file, "/model/elements", rows.size(), fileType.get(), memoryType.get(), + rows.data()); + writeStringAttribute(dataset.get(), "formulation", "FESA-MITC4"); +} + +void writeShellMaterials(const hid_t file, const Domain& domain) { + std::vector rows; + rows.reserve(domain.materials().size()); + for (std::size_t index = 0U; index < domain.materials().size(); ++index) { + const auto& material = domain.materials()[index]; + rows.push_back({ + static_cast(index), + material.name.c_str(), + material.youngsModulus, + material.poissonRatio}); + } + auto stringType = makeUtf8StringType(); + Hdf5Handle fileType{ + requireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(ShellMaterialWriteRow)), + "Unable to create shell material file type."), + H5Tclose}; + Hdf5Handle memoryType{ + requireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(ShellMaterialWriteRow)), + "Unable to create shell material memory type."), + H5Tclose}; + const auto insertFields = [&](const hid_t type, + const hid_t integerType, + const hid_t doubleType) { + requireHdf5( + H5Tinsert(type, "internal_material_id", + HOFFSET(ShellMaterialWriteRow, internalMaterialId), integerType), + "Unable to define shell material ID field."); + requireHdf5( + H5Tinsert(type, "name", + HOFFSET(ShellMaterialWriteRow, name), stringType.get()), + "Unable to define shell material name field."); + requireHdf5( + H5Tinsert(type, "E", HOFFSET(ShellMaterialWriteRow, youngsModulus), + doubleType), + "Unable to define shell material E field."); + requireHdf5( + H5Tinsert(type, "nu", HOFFSET(ShellMaterialWriteRow, poissonRatio), + doubleType), + "Unable to define shell material nu field."); + }; + insertFields(fileType.get(), H5T_STD_U64LE, H5T_IEEE_F64LE); + insertFields(memoryType.get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE); + auto dataset = writeCompoundDataset( + file, "/model/shell/materials", rows.size(), fileType.get(), + memoryType.get(), rows.data()); + writeStringAttribute(dataset.get(), "component_unit_dimensions", "force/length^2,1"); +} + +void writeShellSections(const hid_t file, const Domain& domain) { + std::vector sourceFiles; + sourceFiles.reserve(domain.shellSections().size()); + for (const auto& section : domain.shellSections()) { + sourceFiles.push_back(normalizedPathString(section.location.file)); + } + std::vector rows; + rows.reserve(domain.shellSections().size()); + for (std::size_t index = 0U; index < domain.shellSections().size(); ++index) { + const auto& section = domain.shellSections()[index]; + rows.push_back({ + static_cast(index), + sourceFiles[index].c_str(), + static_cast(section.location.line), + section.name.c_str(), + static_cast(section.materialIndex), + section.thickness}); + } + auto stringType = makeUtf8StringType(); + Hdf5Handle fileType{ + requireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(ShellSectionWriteRow)), + "Unable to create shell section file type."), + H5Tclose}; + Hdf5Handle memoryType{ + requireHdf5Id( + H5Tcreate(H5T_COMPOUND, sizeof(ShellSectionWriteRow)), + "Unable to create shell section memory type."), + H5Tclose}; + const auto insertFields = [&](const hid_t type, + const hid_t integerType, + const hid_t doubleType) { + requireHdf5( + H5Tinsert(type, "internal_section_id", + HOFFSET(ShellSectionWriteRow, internalSectionId), integerType), + "Unable to define shell section ID field."); + requireHdf5( + H5Tinsert(type, "source_file", + HOFFSET(ShellSectionWriteRow, sourceFile), stringType.get()), + "Unable to define shell section source-file field."); + requireHdf5( + H5Tinsert(type, "source_line", + HOFFSET(ShellSectionWriteRow, sourceLine), integerType), + "Unable to define shell section source-line field."); + requireHdf5( + H5Tinsert(type, "source_elset", + HOFFSET(ShellSectionWriteRow, sourceElset), stringType.get()), + "Unable to define shell section ELSET field."); + requireHdf5( + H5Tinsert(type, "material_internal_id", + HOFFSET(ShellSectionWriteRow, materialInternalId), integerType), + "Unable to define shell section material ID field."); + requireHdf5( + H5Tinsert(type, "thickness", + HOFFSET(ShellSectionWriteRow, thickness), doubleType), + "Unable to define shell section thickness field."); + }; + insertFields(fileType.get(), H5T_STD_U64LE, H5T_IEEE_F64LE); + insertFields(memoryType.get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE); + auto dataset = writeCompoundDataset( + file, "/model/shell/sections", rows.size(), fileType.get(), + memoryType.get(), rows.data()); + writeStringAttribute(dataset.get(), "thickness_unit_dimension", "length"); + writeStringAttribute(dataset.get(), "layering", "centered-single-layer"); +} + +void writeShellModelData( + const hid_t file, const Domain& domain, const WriterModelData& modelData) { + std::vector directors; + directors.reserve(domain.nodes().size() * 3U); + std::vector frames; + frames.reserve(domain.nodes().size() * 9U); + for (const auto& frame : domain.shellNodeInitialFrames()) { + directors.insert( + directors.end(), frame.director.begin(), frame.director.end()); + frames.insert(frames.end(), frame.tangentA.begin(), frame.tangentA.end()); + frames.insert(frames.end(), frame.tangentB.begin(), frame.tangentB.end()); + frames.insert(frames.end(), frame.director.begin(), frame.director.end()); + } + writeModelDoubleDataset( + file, "/model/shell/nodal_director", + {static_cast(domain.nodes().size()), 3U}, + directors.data(), directors.size(), "D1,D2,D3", "1,1,1", + "global-cartesian", "nodal"); + writeModelDoubleDataset( + file, "/model/shell/nodal_frame", + {static_cast(domain.nodes().size()), 3U, 3U}, + frames.data(), frames.size(), "X,Y,Z", "1,1,1", + "global-cartesian", "nodal-frame"); + { + Hdf5Handle dataset{ + requireHdf5Id( + H5Dopen2(file, "/model/shell/nodal_frame", H5P_DEFAULT), + "Unable to reopen the nodal-frame dataset."), + H5Dclose}; + writeStringAttribute(dataset.get(), "axis_names", "A,B,D"); + } + + writeShellMaterials(file, domain); + writeShellSections(file, domain); + const std::vector nodalDimensions{ + static_cast(domain.nodes().size()), kDofsPerNode}; + writeUint8Dataset( + file, "/model/nodal_constraint_mask", nodalDimensions, + modelData.constraintMask.data(), modelData.constraintMask.size()); + writeModelDoubleDataset( + file, "/model/prescribed_displacement", nodalDimensions, + modelData.prescribedDisplacement.data(), + modelData.prescribedDisplacement.size(), + "UX,UY,UZ,URX,URY,URZ", + "length,length,length,radian,radian,radian", + "global-cartesian", "nodal-prescribed-value"); + + const double gauss = 1.0 / std::sqrt(3.0); + const std::array locations{ + -gauss, -gauss, + gauss, -gauss, + gauss, gauss, + -gauss, gauss}; + writeModelDoubleDataset( + file, "/model/shell/midsurface_locations", {4U, 2U}, + locations.data(), locations.size(), "XI,ETA", "1,1", + "shell-natural", "midsurface-location"); + const std::array sectionPositions{-1.0, 0.0, 1.0}; + writeModelDoubleDataset( + file, "/model/shell/section_positions", {3U}, + sectionPositions.data(), sectionPositions.size(), "ZETA", "1", + "shell-natural", "section-position"); + { + Hdf5Handle dataset{ + requireHdf5Id( + H5Dopen2(file, "/model/shell/section_positions", H5P_DEFAULT), + "Unable to reopen shell section positions."), + H5Dclose}; + writeStringAttribute(dataset.get(), "position_names", "BOTTOM,MIDDLE,TOP"); + } +} + std::vector flattenEndpointValues( const std::vector& rows, const bool sectionResultants) { @@ -912,6 +1508,155 @@ void writeStress(const hid_t file, const AnalysisState& state) { dataset.get(), "S11", "force/length^2", "beam-local", "section-point"); } +void writeShellResultIdentity( + const hid_t file, + const std::string& path, + const bool usesSectionPositions) { + Hdf5Handle dataset{ + requireHdf5Id( + H5Dopen2(file, path.c_str(), H5P_DEFAULT), + "Unable to reopen a shell result dataset."), + H5Dclose}; + writeStringAttribute( + dataset.get(), + "source_element_type_dataset", + "/model/elements.source_element_type"); + writeStringAttribute(dataset.get(), "internal_formulation", "FESA-MITC4"); + writeStringAttribute( + dataset.get(), + "midsurface_location_dataset", + "/model/shell/midsurface_locations"); + writeStringAttribute( + dataset.get(), + "local_frame_dataset", + "/steps/Step-1/frames/0/element/shell/local_frame"); + if (usesSectionPositions) { + writeStringAttribute( + dataset.get(), + "section_position_dataset", + "/model/shell/section_positions"); + } +} + +void writeShellResultDatasets( + const hid_t file, const Domain& domain, const AnalysisState& state) { + std::vector localFrames; + std::vector generalizedStrains; + std::vector sectionResultants; + std::vector stresses; + localFrames.reserve(state.shellResults().size() * 9U); + generalizedStrains.reserve( + state.shellResults().size() * kShellGeneralizedComponentCount); + sectionResultants.reserve( + state.shellResults().size() * kShellGeneralizedComponentCount); + stresses.reserve( + state.shellResults().size() * kShellSectionPositionCount * + kShellStressComponentCount); + for (const auto& row : state.shellResults()) { + for (const auto& axis : row.localFrame) { + localFrames.insert(localFrames.end(), axis.begin(), axis.end()); + } + generalizedStrains.insert( + generalizedStrains.end(), + row.generalizedStrain.begin(), + row.generalizedStrain.end()); + sectionResultants.insert( + sectionResultants.end(), + row.sectionResultant.begin(), + row.sectionResultant.end()); + for (const auto& position : row.stress) { + stresses.insert( + stresses.end(), + position.components.begin(), + position.components.end()); + } + } + + const hsize_t elementCount = + static_cast(domain.shellElements().size()); + const std::string root = std::string{kStepRoot} + "/element/shell"; + const std::string framePath = root + "/local_frame"; + writeDoubleDataset( + file, framePath, {elementCount, 4U, 3U, 3U}, + localFrames.data(), localFrames.size(), "X,Y,Z", "1,1,1", + "global-cartesian", "shell-local-frame"); + { + Hdf5Handle dataset{ + requireHdf5Id( + H5Dopen2(file, framePath.c_str(), H5P_DEFAULT), + "Unable to reopen shell local-frame results."), + H5Dclose}; + writeStringAttribute(dataset.get(), "axis_names", "E1,E2,E3"); + writeStringAttribute( + dataset.get(), + "source_element_type_dataset", + "/model/elements.source_element_type"); + writeStringAttribute(dataset.get(), "internal_formulation", "FESA-MITC4"); + writeStringAttribute( + dataset.get(), + "midsurface_location_dataset", + "/model/shell/midsurface_locations"); + } + + const std::string strainPath = root + "/generalized_strain"; + writeDoubleDataset( + file, strainPath, {elementCount, 4U, 8U}, + generalizedStrains.data(), generalizedStrains.size(), + "E11,E22,G12,K11,K22,K12,G13,G23", + "1,1,1,1/length,1/length,1/length,1,1", + "shell-local", "midsurface"); + writeShellResultIdentity(file, strainPath, false); + + const std::string resultantPath = root + "/section_resultant"; + writeDoubleDataset( + file, resultantPath, {elementCount, 4U, 8U}, + sectionResultants.data(), sectionResultants.size(), + "N11,N22,N12,M11,M22,M12,Q13,Q23", + "force/length,force/length,force/length,force,force,force,force/length,force/length", + "shell-local", "midsurface"); + writeShellResultIdentity(file, resultantPath, false); + + const std::string stressPath = root + "/stress"; + writeDoubleDataset( + file, stressPath, {elementCount, 4U, 3U, 3U}, + stresses.data(), stresses.size(), "S11,S22,S12", + "force/length^2,force/length^2,force/length^2", + "shell-local", "section-position"); + writeShellResultIdentity(file, stressPath, true); + + const double energy = state.physicalStrainEnergy(); + writeDoubleDataset( + file, std::string{kStepRoot} + "/global/energy", {1U}, + &energy, 1U, "PHYSICAL_STRAIN_ENERGY", "force*length", + "global", "global"); + writeDoubleDataset( + file, std::string{kStepRoot} + "/global/equilibrium", {6U}, + state.equilibrium().data(), state.equilibrium().size(), + "FORCE_1,FORCE_2,FORCE_3,MOMENT_1,MOMENT_2,MOMENT_3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "global-origin"); + const std::string metricsPath = + std::string{kStepRoot} + "/global/verification_metrics"; + writeDoubleDataset( + file, metricsPath, {3U}, state.verificationMetrics().data(), + state.verificationMetrics().size(), + "FREE_RESIDUAL_NORMALIZED,FORCE_BALANCE_NORMALIZED,MOMENT_BALANCE_NORMALIZED", + "1,1,1", "global", "verification"); + { + Hdf5Handle dataset{ + requireHdf5Id( + H5Dopen2(file, metricsPath.c_str(), H5P_DEFAULT), + "Unable to reopen shell verification metrics."), + H5Dclose}; + writeStringAttribute( + dataset.get(), + "metric_definition_ids", + "free-residual-l2-over-max-free-force-l2,force-balance-l2-over-max-force-sum,moment-balance-l2-over-max-moment-sum"); + writeStringAttribute( + dataset.get(), "acceptance_thresholds", "1e-10,1e-10,1e-10"); + } +} + void writeDiagnostics( const hid_t file, const std::vector& inputDiagnostics) { std::vector diagnostics = inputDiagnostics; @@ -1007,6 +1752,11 @@ void writeResultDatasets( "global-cartesian", "nodal"); + if (isShellDomain(domain)) { + writeShellResultDatasets(file, domain, state); + return; + } + const std::vector endpointActionDimensions = { static_cast(domain.elements().size()), kEndpointCount, @@ -1070,7 +1820,7 @@ void writeFile( const Domain& domain, const AnalysisState& state, const std::vector& diagnostics, - const std::vector& localAxes) { + const WriterModelData& modelData) { Hdf5Handle file{ requireHdf5Id( H5Fcreate( @@ -1082,8 +1832,18 @@ void writeFile( (void)createGroup(file.get(), "/model"); (void)createGroup(file.get(), "/steps/Step-1/frames/0/nodal"); (void)createGroup(file.get(), "/steps/Step-1/frames/0/element"); + if (isShellDomain(domain)) { + (void)createGroup(file.get(), "/model/shell"); + (void)createGroup(file.get(), "/steps/Step-1/frames/0/element/shell"); + (void)createGroup(file.get(), "/steps/Step-1/frames/0/global"); + } writeNodes(file.get(), domain); - writeElements(file.get(), domain, localAxes); + if (isShellDomain(domain)) { + writeShellElements(file.get(), domain); + writeShellModelData(file.get(), domain, modelData); + } else { + writeBeamElements(file.get(), domain, modelData.beamLocalAxes); + } writeResultDatasets(file.get(), domain, state); writeDiagnostics(file.get(), diagnostics); requireHdf5( @@ -1185,11 +1945,18 @@ void requirePortableCompoundMember( H5Tclose}; const bool isUint64 = name == "internal_node_id" || name == "internal_element_id" || + name == "internal_material_id" || name == "internal_section_id" || + name == "shell_section_internal_id" || + name == "material_internal_id" || name == "source_line" || name == "gauss_point_index" || name == "section_point_index" || name == "line"; - const bool isFloat64 = name == "x1" || name == "x2" || name == "S11"; + const bool isFloat64 = + name == "x1" || name == "x2" || name == "S11" || + name == "E" || name == "nu" || name == "thickness"; const bool isString = name == "instance_name" || name == "source_label" || + name == "source_element_type" || name == "internal_formulation" || + name == "name" || name == "source_file" || name == "source_elset" || name == "source" || name == "severity" || name == "code" || name == "file" || name == "keyword" || name == "entity_identity" || name == "message"; @@ -1236,9 +2003,10 @@ void requirePortableCompoundMember( "Unable to inspect a compound HDF5 array base type."), H5Tclose}; if (name == "node_internal_ids") { - if (dimensions != std::vector{2U} || + if ((dimensions != std::vector{2U} && + dimensions != std::vector{4U}) || H5Tequal(baseType.get(), H5T_STD_U64LE) <= 0) { - throw Hdf5Failure{"Element connectivity is not uint64[2]."}; + throw Hdf5Failure{"Element connectivity is not uint64[2] or uint64[4]."}; } } else if (name == "coordinates") { if (dimensions != std::vector{3U} || @@ -1270,6 +2038,31 @@ void requireResultAttributes( requireUint64Attribute(dataset, "frame_index", 0U); } +void requireShellResultIdentity( + const hid_t file, const char* path, const bool usesSectionPositions) { + auto dataset = openDatasetForCheck(file, path); + requireStringAttribute( + dataset.get(), + "source_element_type_dataset", + "/model/elements.source_element_type"); + requireStringAttribute( + dataset.get(), "internal_formulation", "FESA-MITC4"); + requireStringAttribute( + dataset.get(), + "midsurface_location_dataset", + "/model/shell/midsurface_locations"); + requireStringAttribute( + dataset.get(), + "local_frame_dataset", + "/steps/Step-1/frames/0/element/shell/local_frame"); + if (usesSectionPositions) { + requireStringAttribute( + dataset.get(), + "section_position_dataset", + "/model/shell/section_positions"); + } +} + void requireDoubleDataset( const hid_t file, const char* path, @@ -1293,6 +2086,116 @@ void requireDoubleDataset( } requireResultAttributes( dataset.get(), componentNames, componentUnits, coordinateSystem, location); + std::size_t valueCount = 1U; + for (const hsize_t dimension : expectedDimensions) { + std::size_t next = 0U; + if (!sizeProductFits( + valueCount, static_cast(dimension), next)) { + throw Hdf5Failure{"A result dataset shape overflows size_t."}; + } + valueCount = next; + } + std::vector values(valueCount); + if (!values.empty()) { + requireHdf5( + H5Dread( + dataset.get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), + "Unable to read a result dataset during self-check."); + } + if (!std::all_of(values.begin(), values.end(), [](const double value) { + return std::isfinite(value); + })) { + throw Hdf5Failure{"A result dataset contains a nonfinite value."}; + } +} + +void requireModelDoubleDataset( + const hid_t file, + const char* path, + const std::vector& expectedDimensions, + const char* componentNames, + const char* componentUnits, + const char* coordinateSystem, + const char* location, + const std::vector* expectedValues = nullptr) { + auto dataset = openDatasetForCheck(file, path); + if (checkedDimensions(dataset.get()) != expectedDimensions) { + throw Hdf5Failure{"A model floating-point dataset has the wrong shape."}; + } + Hdf5Handle type{ + requireHdf5Id( + H5Dget_type(dataset.get()), "Unable to inspect a model dataset type."), + H5Tclose}; + if (H5Tget_class(type.get()) != H5T_FLOAT || + H5Tget_size(type.get()) != sizeof(double) || + H5Tequal(type.get(), H5T_IEEE_F64LE) <= 0) { + throw Hdf5Failure{"A model dataset is not IEEE-754 float64 little-endian."}; + } + requireStringAttribute(dataset.get(), "component_names", componentNames); + requireStringAttribute( + dataset.get(), "component_unit_dimensions", componentUnits); + requireStringAttribute(dataset.get(), "coordinate_system", coordinateSystem); + requireStringAttribute(dataset.get(), "location", location); + std::size_t valueCount = 1U; + for (const hsize_t dimension : expectedDimensions) { + std::size_t next = 0U; + if (!sizeProductFits( + valueCount, static_cast(dimension), next)) { + throw Hdf5Failure{"A model dataset shape overflows size_t."}; + } + valueCount = next; + } + std::vector values(valueCount); + if (!values.empty()) { + requireHdf5( + H5Dread( + dataset.get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), + "Unable to read a model dataset during self-check."); + } + if (!std::all_of(values.begin(), values.end(), [](const double value) { + return std::isfinite(value); + })) { + throw Hdf5Failure{"A model dataset contains a nonfinite value."}; + } + if (expectedValues != nullptr && values != *expectedValues) { + throw Hdf5Failure{"A fixed model dataset has the wrong value order."}; + } +} + +void requireUint8Dataset( + const hid_t file, + const char* path, + const std::vector& expectedDimensions, + const std::vector& expectedValues) { + auto dataset = openDatasetForCheck(file, path); + if (checkedDimensions(dataset.get()) != expectedDimensions) { + throw Hdf5Failure{"A uint8 HDF5 dataset has the wrong shape."}; + } + Hdf5Handle type{ + requireHdf5Id( + H5Dget_type(dataset.get()), "Unable to inspect a uint8 dataset type."), + H5Tclose}; + if (H5Tget_class(type.get()) != H5T_INTEGER || + H5Tget_size(type.get()) != sizeof(std::uint8_t) || + H5Tget_sign(type.get()) != H5T_SGN_NONE || + H5Tequal(type.get(), H5T_STD_U8LE) <= 0) { + throw Hdf5Failure{"A constraint mask is not portable uint8."}; + } + std::vector values(expectedValues.size()); + if (!values.empty()) { + requireHdf5( + H5Dread( + dataset.get(), H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), + "Unable to read the constraint mask during self-check."); + } + if (values != expectedValues) { + throw Hdf5Failure{"The constraint mask has the wrong value order."}; + } + requireStringAttribute(dataset.get(), "component_names", "UX,UY,UZ,URX,URY,URZ"); + requireStringAttribute(dataset.get(), "value_meaning", "0=free,1=constrained"); } void requireCompoundDataset( @@ -1333,7 +2236,8 @@ void selfCheckFile( const std::filesystem::path& path, const Domain& domain, const AnalysisState& state, - const std::size_t diagnosticCount) { + const std::size_t diagnosticCount, + const WriterModelData& modelData) { if (H5Fis_hdf5(path.string().c_str()) <= 0) { throw Hdf5Failure{"The temporary output is not an HDF5 file."}; } @@ -1342,27 +2246,41 @@ void selfCheckFile( H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), "Unable to reopen the temporary HDF5 file read-only."), H5Fclose}; - { Hdf5Handle metadata{ requireHdf5Id( H5Gopen2(file.get(), "/metadata", H5P_DEFAULT), "The HDF5 metadata group is missing."), H5Gclose}; requireUint64Attribute(metadata.get(), "schema_version", 0U); - requireStringAttribute( - metadata.get(), "feature_id", "linear-static-3d-euler-beam"); requireStringAttribute( metadata.get(), "solver_version", std::string{solverVersion()}); requireStringAttribute( metadata.get(), "source_input_identity", sourceInputIdentity(domain)); requireStringAttribute( metadata.get(), "unit_system_label", "user-consistent-unspecified"); - requireStringAttribute( - metadata.get(), - "coordinate_convention", - "global-cartesian; beam-local=(t,n1,t-cross-n1)"); - requireStringAttribute( - metadata.get(), "element_formulation", "B33-3D-Euler-Bernoulli"); + if (isShellDomain(domain)) { + requireStringAttribute( + metadata.get(), "feature_id", "linear-static-mitc4-shell"); + requireStringAttribute( + metadata.get(), + "coordinate_convention", + "global-cartesian; shell-local=(e1,e2,e3); positive-thickness=+zeta"); + requireStringAttribute( + metadata.get(), "internal_formulation", "FESA-MITC4"); + requireStringAttribute( + metadata.get(), + "integration_rule", + "2x2x2-gauss; mitc4-edge-midpoint-shear"); + } else { + requireStringAttribute( + metadata.get(), "feature_id", "linear-static-3d-euler-beam"); + requireStringAttribute( + metadata.get(), + "coordinate_convention", + "global-cartesian; beam-local=(t,n1,t-cross-n1)"); + requireStringAttribute( + metadata.get(), "element_formulation", "B33-3D-Euler-Bernoulli"); + } requireStringAttribute(metadata.get(), "step_name", kStepName); requireUint64Attribute(metadata.get(), "frame_index", 0U); @@ -1374,26 +2292,8 @@ void selfCheckFile( auto nodes = openDatasetForCheck(file.get(), "/model/nodes"); requireStringAttribute(nodes.get(), "coordinate_system", "global-cartesian"); requireStringAttribute(nodes.get(), "units_label", "length"); - requireCompoundDataset( - file.get(), - "/model/elements", - static_cast(domain.elements().size()), - {"internal_element_id", "instance_name", "source_label", - "node_internal_ids", "local_axes"}); - auto elements = openDatasetForCheck(file.get(), "/model/elements"); - requireStringAttribute( - elements.get(), "formulation", "B33-3D-Euler-Bernoulli"); - const std::vector nodalDimensions = { static_cast(domain.nodes().size()), kDofsPerNode}; - const std::vector endDimensions = { - static_cast(domain.elements().size()), - kEndpointCount, - kEndActionComponentCount}; - const std::vector generalizedDimensions = { - static_cast(domain.elements().size()), - kGaussPointCount, - kGeneralizedComponentCount}; requireDoubleDataset( file.get(), "/steps/Step-1/frames/0/nodal/displacement", @@ -1410,55 +2310,204 @@ void selfCheckFile( "force,force,force,force*length,force*length,force*length", "global-cartesian", "nodal"); - requireDoubleDataset( - file.get(), - "/steps/Step-1/frames/0/element/end_force_local", - endDimensions, - "FX,FY,FZ,MX,MY,MZ", - "force,force,force,force*length,force*length,force*length", - "beam-local", - "endpoint-outward-action"); - requireDoubleDataset( - file.get(), - "/steps/Step-1/frames/0/element/section_resultant", - generalizedDimensions, - "N,T,My,Mz", - "force,force*length,force*length,force*length", - "beam-local", - "endpoint-positive-local-x-section-cut"); - requireDoubleDataset( - file.get(), - "/steps/Step-1/frames/0/element/generalized_strain", - generalizedDimensions, - "epsilon0,kappa_x,kappa_y,kappa_z", - "1,1/length,1/length,1/length", - "beam-local", - "integration-point"); - requireDoubleDataset( - file.get(), - "/steps/Step-1/frames/0/element/generalized_resultant", - generalizedDimensions, - "N,T,My,Mz", - "force,force*length,force*length,force*length", - "beam-local", - "integration-point"); - requireCompoundDataset( - file.get(), - "/steps/Step-1/frames/0/element/stress_s11", - static_cast(state.stressResults().size()), - {"internal_element_id", "gauss_point_index", "section_point_index", - "x1", "x2", "source", "S11"}); - auto stress = openDatasetForCheck( - file.get(), "/steps/Step-1/frames/0/element/stress_s11"); - requireResultAttributes( - stress.get(), "S11", "force/length^2", "beam-local", "section-point"); + + if (isShellDomain(domain)) { + requireCompoundDataset( + file.get(), + "/model/elements", + static_cast(domain.shellElements().size()), + {"internal_element_id", "instance_name", "source_label", + "source_element_type", "internal_formulation", "node_internal_ids", + "shell_section_internal_id", "material_internal_id"}); + auto elements = openDatasetForCheck(file.get(), "/model/elements"); + requireStringAttribute(elements.get(), "formulation", "FESA-MITC4"); + requireCompoundDataset( + file.get(), + "/model/shell/materials", + static_cast(domain.materials().size()), + {"internal_material_id", "name", "E", "nu"}); + requireCompoundDataset( + file.get(), + "/model/shell/sections", + static_cast(domain.shellSections().size()), + {"internal_section_id", "source_file", "source_line", "source_elset", + "material_internal_id", "thickness"}); + + std::vector directors; + std::vector frames; + directors.reserve(domain.nodes().size() * 3U); + frames.reserve(domain.nodes().size() * 9U); + for (const auto& frame : domain.shellNodeInitialFrames()) { + directors.insert( + directors.end(), frame.director.begin(), frame.director.end()); + frames.insert(frames.end(), frame.tangentA.begin(), frame.tangentA.end()); + frames.insert(frames.end(), frame.tangentB.begin(), frame.tangentB.end()); + frames.insert(frames.end(), frame.director.begin(), frame.director.end()); + } + requireModelDoubleDataset( + file.get(), "/model/shell/nodal_director", + {static_cast(domain.nodes().size()), 3U}, + "D1,D2,D3", "1,1,1", "global-cartesian", "nodal", &directors); + requireModelDoubleDataset( + file.get(), "/model/shell/nodal_frame", + {static_cast(domain.nodes().size()), 3U, 3U}, + "X,Y,Z", "1,1,1", "global-cartesian", "nodal-frame", &frames); + auto nodalFrame = openDatasetForCheck( + file.get(), "/model/shell/nodal_frame"); + requireStringAttribute(nodalFrame.get(), "axis_names", "A,B,D"); + requireUint8Dataset( + file.get(), "/model/nodal_constraint_mask", nodalDimensions, + modelData.constraintMask); + requireModelDoubleDataset( + file.get(), "/model/prescribed_displacement", nodalDimensions, + "UX,UY,UZ,URX,URY,URZ", + "length,length,length,radian,radian,radian", + "global-cartesian", "nodal-prescribed-value", + &modelData.prescribedDisplacement); + const double gauss = 1.0 / std::sqrt(3.0); + const std::vector locations{ + -gauss, -gauss, gauss, -gauss, + gauss, gauss, -gauss, gauss}; + requireModelDoubleDataset( + file.get(), "/model/shell/midsurface_locations", {4U, 2U}, + "XI,ETA", "1,1", "shell-natural", "midsurface-location", + &locations); + const std::vector sectionPositions{-1.0, 0.0, 1.0}; + requireModelDoubleDataset( + file.get(), "/model/shell/section_positions", {3U}, + "ZETA", "1", "shell-natural", "section-position", + §ionPositions); + auto positions = openDatasetForCheck( + file.get(), "/model/shell/section_positions"); + requireStringAttribute( + positions.get(), "position_names", "BOTTOM,MIDDLE,TOP"); + + const hsize_t elementCount = + static_cast(domain.shellElements().size()); + const std::string shellRoot = std::string{kStepRoot} + "/element/shell"; + requireDoubleDataset( + file.get(), (shellRoot + "/local_frame").c_str(), + {elementCount, 4U, 3U, 3U}, "X,Y,Z", "1,1,1", + "global-cartesian", "shell-local-frame"); + auto localFrame = openDatasetForCheck( + file.get(), (shellRoot + "/local_frame").c_str()); + requireStringAttribute(localFrame.get(), "axis_names", "E1,E2,E3"); + requireStringAttribute( + localFrame.get(), "internal_formulation", "FESA-MITC4"); + requireStringAttribute( + localFrame.get(), + "source_element_type_dataset", + "/model/elements.source_element_type"); + requireStringAttribute( + localFrame.get(), + "midsurface_location_dataset", + "/model/shell/midsurface_locations"); + requireDoubleDataset( + file.get(), (shellRoot + "/generalized_strain").c_str(), + {elementCount, 4U, 8U}, + "E11,E22,G12,K11,K22,K12,G13,G23", + "1,1,1,1/length,1/length,1/length,1,1", + "shell-local", "midsurface"); + requireShellResultIdentity( + file.get(), (shellRoot + "/generalized_strain").c_str(), false); + requireDoubleDataset( + file.get(), (shellRoot + "/section_resultant").c_str(), + {elementCount, 4U, 8U}, + "N11,N22,N12,M11,M22,M12,Q13,Q23", + "force/length,force/length,force/length,force,force,force,force/length,force/length", + "shell-local", "midsurface"); + requireShellResultIdentity( + file.get(), (shellRoot + "/section_resultant").c_str(), false); + requireDoubleDataset( + file.get(), (shellRoot + "/stress").c_str(), + {elementCount, 4U, 3U, 3U}, "S11,S22,S12", + "force/length^2,force/length^2,force/length^2", + "shell-local", "section-position"); + requireShellResultIdentity( + file.get(), (shellRoot + "/stress").c_str(), true); + requireDoubleDataset( + file.get(), "/steps/Step-1/frames/0/global/energy", {1U}, + "PHYSICAL_STRAIN_ENERGY", "force*length", "global", "global"); + requireDoubleDataset( + file.get(), "/steps/Step-1/frames/0/global/equilibrium", {6U}, + "FORCE_1,FORCE_2,FORCE_3,MOMENT_1,MOMENT_2,MOMENT_3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "global-origin"); + requireDoubleDataset( + file.get(), "/steps/Step-1/frames/0/global/verification_metrics", {3U}, + "FREE_RESIDUAL_NORMALIZED,FORCE_BALANCE_NORMALIZED,MOMENT_BALANCE_NORMALIZED", + "1,1,1", "global", "verification"); + auto metrics = openDatasetForCheck( + file.get(), "/steps/Step-1/frames/0/global/verification_metrics"); + requireStringAttribute( + metrics.get(), + "metric_definition_ids", + "free-residual-l2-over-max-free-force-l2,force-balance-l2-over-max-force-sum,moment-balance-l2-over-max-moment-sum"); + requireStringAttribute( + metrics.get(), "acceptance_thresholds", "1e-10,1e-10,1e-10"); + for (const char* forbidden : { + "/steps/Step-1/frames/0/element/shell/drilling", + "/steps/Step-1/frames/0/element/shell/drilling_energy", + "/steps/Step-1/frames/0/element/shell/S33", + "/steps/Step-1/frames/0/element/shell/S13", + "/steps/Step-1/frames/0/element/shell/S23"}) { + if (H5Lexists(file.get(), forbidden, H5P_DEFAULT) != 0) { + throw Hdf5Failure{"A forbidden shell result path exists."}; + } + } + } else { + requireCompoundDataset( + file.get(), + "/model/elements", + static_cast(domain.elements().size()), + {"internal_element_id", "instance_name", "source_label", + "node_internal_ids", "local_axes"}); + auto elements = openDatasetForCheck(file.get(), "/model/elements"); + requireStringAttribute( + elements.get(), "formulation", "B33-3D-Euler-Bernoulli"); + const std::vector endDimensions = { + static_cast(domain.elements().size()), + kEndpointCount, + kEndActionComponentCount}; + const std::vector generalizedDimensions = { + static_cast(domain.elements().size()), + kGaussPointCount, + kGeneralizedComponentCount}; + requireDoubleDataset( + file.get(), "/steps/Step-1/frames/0/element/end_force_local", + endDimensions, "FX,FY,FZ,MX,MY,MZ", + "force,force,force,force*length,force*length,force*length", + "beam-local", "endpoint-outward-action"); + requireDoubleDataset( + file.get(), "/steps/Step-1/frames/0/element/section_resultant", + generalizedDimensions, "N,T,My,Mz", + "force,force*length,force*length,force*length", + "beam-local", "endpoint-positive-local-x-section-cut"); + requireDoubleDataset( + file.get(), "/steps/Step-1/frames/0/element/generalized_strain", + generalizedDimensions, "epsilon0,kappa_x,kappa_y,kappa_z", + "1,1/length,1/length,1/length", "beam-local", "integration-point"); + requireDoubleDataset( + file.get(), "/steps/Step-1/frames/0/element/generalized_resultant", + generalizedDimensions, "N,T,My,Mz", + "force,force*length,force*length,force*length", + "beam-local", "integration-point"); + requireCompoundDataset( + file.get(), "/steps/Step-1/frames/0/element/stress_s11", + static_cast(state.stressResults().size()), + {"internal_element_id", "gauss_point_index", "section_point_index", + "x1", "x2", "source", "S11"}); + auto stress = openDatasetForCheck( + file.get(), "/steps/Step-1/frames/0/element/stress_s11"); + requireResultAttributes( + stress.get(), "S11", "force/length^2", "beam-local", "section-point"); + } requireCompoundDataset( file.get(), "/diagnostics", static_cast(diagnosticCount), {"severity", "code", "file", "line", "keyword", "entity_identity", "message"}); - } requireHdf5( file.closeChecked(), "Unable to close the read-only HDF5 schema self-check handle."); @@ -1519,10 +2568,10 @@ Status Hdf5ResultsWriter::write( const Domain& domain, const AnalysisState& state, const std::vector& diagnostics) { - std::vector localAxes; + WriterModelData modelData; try { const Status validation = validateWriterInput( - outputPath, domain, state, diagnostics, localAxes); + outputPath, domain, state, diagnostics, modelData); if (!validation.isOk()) { return validation; } @@ -1530,11 +2579,11 @@ Status Hdf5ResultsWriter::write( Hdf5ErrorSilencer silenceBackendErrors; const std::filesystem::path temporaryPath = makeTemporaryPath(outputPath); TemporaryFileGuard cleanup{temporaryPath}; - writeFile(temporaryPath, domain, state, diagnostics, localAxes); + writeFile(temporaryPath, domain, state, diagnostics, modelData); // Close occurs when writeFile returns; only the reopened read-only file // can approve this temp artifact for authoritative replacement. selfCheckFile( - temporaryPath, domain, state, diagnostics.size()); + temporaryPath, domain, state, diagnostics.size(), modelData); if (!finalizeFile(temporaryPath, outputPath)) { return outputFailure( "hdf5-finalization-failure", diff --git a/tests/unit/io/hdf5/hdf5_results_writer_test.cpp b/tests/unit/io/hdf5/hdf5_results_writer_test.cpp index 35a7f9a..f15f5d4 100644 --- a/tests/unit/io/hdf5/hdf5_results_writer_test.cpp +++ b/tests/unit/io/hdf5/hdf5_results_writer_test.cpp @@ -215,6 +215,122 @@ WriterFixture makeFixture( return {std::move(domain), std::move(dofs), std::move(state)}; } +fesa::ModelDefinition makeShellDefinition(const std::filesystem::path& source) { + fesa::ModelDefinition definition{}; + definition.sourcePath = source; + definition.sourceContentIdentity = "fnv1a64:fedcba9876543210"; + definition.nodes = { + {{"Shell-1", 11, "11"}, {-1.0, -1.0, 0.0}, {source, 10U}}, + {{"Shell-1", 12, "12"}, {1.0, -1.0, 0.0}, {source, 11U}}, + {{"Shell-1", 13, "13"}, {1.0, 1.0, 0.0}, {source, 12U}}, + {{"Shell-1", 14, "14"}, {-1.0, 1.0, 0.0}, {source, 13U}}}; + definition.materials = { + {"ShellSteel", 210.0e9, 0.3, {source, 20U}}}; + definition.shellSections = { + {"PlateSet", 0.02, 0U, {source, 30U}}}; + definition.shellElements = {{ + {"Shell-1", 401, "401"}, + fesa::ShellSourceElementType::s4r, + {0U, 1U, 2U, 3U}, + 0U, + 0U, + {source, 40U}}}; + for (std::size_t node = 0U; node < definition.nodes.size(); ++node) { + definition.shellNodeInitialFrames.push_back({ + static_cast(node), + {0.0, 0.0, 1.0}, + {1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}}); + } + definition.nodeSets = { + {"Fixed", {}, {0U}, {source, 50U}}}; + definition.steps = {{ + "Step-1", + {{"Fixed", 1, 3, 0.0, {source, 60U}}, + {"Fixed", 4, 4, 0.125, {source, 61U}}}, + {}, + 0.1, + 1.0, + 0.01, + 1.0, + {source, 59U}}}; + return definition; +} + +WriterFixture makeShellFixture(const std::filesystem::path& source) { + auto domainResult = fesa::Domain::create(makeShellDefinition(source)); + if (!domainResult.hasValue()) { + throw std::runtime_error{"Shell writer fixture Domain construction failed."}; + } + auto domain = std::make_unique( + std::move(domainResult.value())); + auto modelResult = fesa::AnalysisModel::create(*domain); + if (!modelResult.hasValue()) { + throw std::runtime_error{"Shell writer fixture AnalysisModel construction failed."}; + } + const fesa::AnalysisModel model = std::move(modelResult.value()); + auto dofsResult = fesa::DofManager::create(model); + if (!dofsResult.hasValue()) { + throw std::runtime_error{"Shell writer fixture DofManager construction failed."}; + } + auto dofs = std::make_unique( + std::move(dofsResult.value())); + auto state = std::make_unique( + fesa::AnalysisState::create(*dofs, {"Step-1", 0U})); + for (std::size_t index = 0U; index < state->displacement().size(); ++index) { + state->displacement()[index] = 0.01 * static_cast(index + 1U); + state->externalForce()[index] = 10.0 + static_cast(index); + state->internalForce()[index] = 20.0 + static_cast(index); + state->residual()[index] = 30.0 + static_cast(index); + state->reaction()[index] = 40.0 + static_cast(index); + } + + const double gauss = 1.0 / std::sqrt(3.0); + const std::array, 4> coordinates{{ + {-gauss, -gauss}, + {gauss, -gauss}, + {gauss, gauss}, + {-gauss, gauss}}}; + const std::array locations{ + fesa::ShellMidsurfaceLocation::gp1, + fesa::ShellMidsurfaceLocation::gp2, + fesa::ShellMidsurfaceLocation::gp3, + fesa::ShellMidsurfaceLocation::gp4}; + fesa::ShellStateCandidate candidate{}; + for (std::size_t point = 0U; point < locations.size(); ++point) { + const double base = 100.0 * static_cast(point + 1U); + candidate.rows.push_back({ + 0U, + locations[point], + coordinates[point], + {{{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}}}, + {base + 1.0, base + 2.0, base + 3.0, base + 4.0, + base + 5.0, base + 6.0, base + 7.0, base + 8.0}, + {base + 11.0, base + 12.0, base + 13.0, base + 14.0, + base + 15.0, base + 16.0, base + 17.0, base + 18.0}, + {{{fesa::ShellSectionPosition::bottom, + -1.0, + {base + 21.0, base + 22.0, base + 23.0}}, + {fesa::ShellSectionPosition::middle, + 0.0, + {base + 24.0, base + 25.0, base + 26.0}}, + {fesa::ShellSectionPosition::top, + 1.0, + {base + 27.0, base + 28.0, base + 29.0}}}}}); + } + candidate.physicalStrainEnergy = 123.5; + candidate.equilibrium = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + candidate.verificationMetrics = {1.0e-13, 2.0e-13, 3.0e-13}; + const fesa::Status commit = state->commitShellResults( + {0U}, std::move(candidate)); + if (!commit.isOk()) { + throw std::runtime_error{"Shell writer fixture state commit failed."}; + } + return {std::move(domain), std::move(dofs), std::move(state)}; +} + Hdf5Handle openFile(const std::filesystem::path& path) { const hid_t file = H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); if (file < 0) { @@ -267,6 +383,30 @@ std::vector readDoubleDataset( return values; } +std::vector readUint8Dataset( + const hid_t file, const std::string& path) { + const auto dimensions = datasetDimensions(file, path); + std::size_t valueCount = 1U; + for (const hsize_t dimension : dimensions) { + valueCount *= static_cast(dimension); + } + const auto dataset = openDataset(file, path); + Hdf5Handle type{H5Dget_type(dataset.get()), H5Tclose}; + if (type.get() < 0 || H5Tget_class(type.get()) != H5T_INTEGER || + H5Tget_size(type.get()) != sizeof(std::uint8_t) || + H5Tget_sign(type.get()) != H5T_SGN_NONE || + H5Tequal(type.get(), H5T_STD_U8LE) <= 0) { + throw std::runtime_error{"Expected a portable uint8 HDF5 dataset."}; + } + std::vector values(valueCount); + if (!values.empty() && + H5Dread(dataset.get(), H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()) < 0) { + throw std::runtime_error{"Unable to read uint8 HDF5 dataset."}; + } + return values; +} + std::string readStringAttribute(const hid_t object, const char* name) { Hdf5Handle attribute{H5Aopen(object, name, H5P_DEFAULT), H5Aclose}; Hdf5Handle type{H5Aget_type(attribute.get()), H5Tclose}; @@ -380,6 +520,25 @@ void expectCompoundMembers( } } +void expectCompoundMemberNames( + const hid_t file, + const std::string& path, + const std::vector& expectedNames) { + const auto dataset = openDataset(file, path); + Hdf5Handle type{H5Dget_type(dataset.get()), H5Tclose}; + ASSERT_EQ(H5Tget_class(type.get()), H5T_COMPOUND); + ASSERT_EQ( + H5Tget_nmembers(type.get()), static_cast(expectedNames.size())); + for (std::size_t index = 0U; index < expectedNames.size(); ++index) { + char* rawName = H5Tget_member_name( + type.get(), static_cast(index)); + ASSERT_NE(rawName, nullptr); + const std::string actualName{rawName}; + (void)H5free_memory(rawName); + EXPECT_EQ(actualName, expectedNames[index]); + } +} + void expectNumericDataset( const hid_t file, const std::string& path, @@ -1003,3 +1162,213 @@ TEST(Hdf5ResultsWriter, SuccessfullyReplacesExistingFinal) { ASSERT_GE(metadata.get(), 0); EXPECT_EQ(readUint64Attribute(metadata.get(), "schema_version"), 0U); } + +// MITC4-H5-001 +TEST(Hdf5ResultsWriter, WritesExactShellMetadataAndModelIdentity) { + TempDirectory directory{"shell-model"}; + const auto source = directory.path() / "shell.inp"; + auto fixture = makeShellFixture(source); + const auto output = directory.path() / "results.h5"; + + fesa::Hdf5ResultsWriter writer; + ASSERT_TRUE(writer.write(output, *fixture.domain, *fixture.state, {}).isOk()); + const auto file = openFile(output); + + Hdf5Handle metadata{ + H5Gopen2(file.get(), "/metadata", H5P_DEFAULT), H5Gclose}; + ASSERT_GE(metadata.get(), 0); + EXPECT_EQ(readUint64Attribute(metadata.get(), "schema_version"), 0U); + EXPECT_EQ( + readStringAttribute(metadata.get(), "feature_id"), + "linear-static-mitc4-shell"); + EXPECT_EQ( + readStringAttribute(metadata.get(), "coordinate_convention"), + "global-cartesian; shell-local=(e1,e2,e3); positive-thickness=+zeta"); + EXPECT_EQ( + readStringAttribute(metadata.get(), "internal_formulation"), + "FESA-MITC4"); + EXPECT_EQ( + readStringAttribute(metadata.get(), "integration_rule"), + "2x2x2-gauss; mitc4-edge-midpoint-shear"); + + EXPECT_EQ( + datasetDimensions(file.get(), "/model/elements"), + std::vector({1U})); + expectCompoundMemberNames( + file.get(), + "/model/elements", + {"internal_element_id", "instance_name", "source_label", + "source_element_type", "internal_formulation", "node_internal_ids", + "shell_section_internal_id", "material_internal_id"}); + const auto elements = openDataset(file.get(), "/model/elements"); + EXPECT_EQ( + readStringAttribute(elements.get(), "formulation"), "FESA-MITC4"); + + EXPECT_EQ( + datasetDimensions(file.get(), "/model/shell/nodal_director"), + std::vector({4U, 3U})); + EXPECT_EQ( + readDoubleDataset(file.get(), "/model/shell/nodal_director"), + std::vector({0.0, 0.0, 1.0, + 0.0, 0.0, 1.0, + 0.0, 0.0, 1.0, + 0.0, 0.0, 1.0})); + EXPECT_EQ( + datasetDimensions(file.get(), "/model/shell/nodal_frame"), + std::vector({4U, 3U, 3U})); + expectCompoundMemberNames( + file.get(), "/model/shell/materials", + {"internal_material_id", "name", "E", "nu"}); + expectCompoundMemberNames( + file.get(), "/model/shell/sections", + {"internal_section_id", "source_file", "source_line", "source_elset", + "material_internal_id", "thickness"}); + + EXPECT_EQ( + datasetDimensions(file.get(), "/model/nodal_constraint_mask"), + std::vector({4U, 6U})); + const auto mask = readUint8Dataset(file.get(), "/model/nodal_constraint_mask"); + ASSERT_EQ(mask.size(), 24U); + EXPECT_EQ(mask[0U], 1U); + EXPECT_EQ(mask[1U], 1U); + EXPECT_EQ(mask[2U], 1U); + EXPECT_EQ(mask[3U], 1U); + EXPECT_EQ(mask[4U], 0U); + EXPECT_EQ(mask[5U], 0U); + const auto prescribed = readDoubleDataset( + file.get(), "/model/prescribed_displacement"); + ASSERT_EQ(prescribed.size(), 24U); + EXPECT_DOUBLE_EQ(prescribed[0U], 0.0); + EXPECT_DOUBLE_EQ(prescribed[3U], 0.125); + EXPECT_DOUBLE_EQ(prescribed[4U], 0.0); + EXPECT_EQ( + readDoubleDataset(file.get(), "/model/shell/section_positions"), + std::vector({-1.0, 0.0, 1.0})); + const auto locations = readDoubleDataset( + file.get(), "/model/shell/midsurface_locations"); + ASSERT_EQ(locations.size(), 8U); + const double gauss = 1.0 / std::sqrt(3.0); + EXPECT_DOUBLE_EQ(locations[0U], -gauss); + EXPECT_DOUBLE_EQ(locations[1U], -gauss); + EXPECT_DOUBLE_EQ(locations[6U], -gauss); + EXPECT_DOUBLE_EQ(locations[7U], gauss); +} + +// MITC4-H5-002 +TEST(Hdf5ResultsWriter, WritesExactMandatoryShellResultInventory) { + TempDirectory directory{"shell-results"}; + auto fixture = makeShellFixture(directory.path() / "shell.inp"); + const auto output = directory.path() / "results.h5"; + + fesa::Hdf5ResultsWriter writer; + ASSERT_TRUE(writer.write(output, *fixture.domain, *fixture.state, {}).isOk()); + const auto file = openFile(output); + const std::string shellRoot = std::string{kStepRoot} + "/element/shell"; + expectNumericDataset( + file.get(), shellRoot + "/local_frame", {1U, 4U, 3U, 3U}, + "X,Y,Z", "1,1,1", "global-cartesian", "shell-local-frame"); + expectNumericDataset( + file.get(), shellRoot + "/generalized_strain", {1U, 4U, 8U}, + "E11,E22,G12,K11,K22,K12,G13,G23", + "1,1,1,1/length,1/length,1/length,1,1", + "shell-local", "midsurface"); + expectNumericDataset( + file.get(), shellRoot + "/section_resultant", {1U, 4U, 8U}, + "N11,N22,N12,M11,M22,M12,Q13,Q23", + "force/length,force/length,force/length,force,force,force,force/length,force/length", + "shell-local", "midsurface"); + expectNumericDataset( + file.get(), shellRoot + "/stress", {1U, 4U, 3U, 3U}, + "S11,S22,S12", + "force/length^2,force/length^2,force/length^2", + "shell-local", "section-position"); + + const auto strain = readDoubleDataset(file.get(), shellRoot + "/generalized_strain"); + ASSERT_EQ(strain.size(), 32U); + EXPECT_DOUBLE_EQ(strain.front(), 101.0); + EXPECT_DOUBLE_EQ(strain.back(), 408.0); + const auto stress = readDoubleDataset(file.get(), shellRoot + "/stress"); + ASSERT_EQ(stress.size(), 36U); + EXPECT_DOUBLE_EQ(stress.front(), 121.0); + EXPECT_DOUBLE_EQ(stress.back(), 429.0); + + expectNumericDataset( + file.get(), std::string{kStepRoot} + "/global/energy", {1U}, + "PHYSICAL_STRAIN_ENERGY", "force*length", "global", "global"); + expectNumericDataset( + file.get(), std::string{kStepRoot} + "/global/equilibrium", {6U}, + "FORCE_1,FORCE_2,FORCE_3,MOMENT_1,MOMENT_2,MOMENT_3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "global-origin"); + expectNumericDataset( + file.get(), std::string{kStepRoot} + "/global/verification_metrics", {3U}, + "FREE_RESIDUAL_NORMALIZED,FORCE_BALANCE_NORMALIZED,MOMENT_BALANCE_NORMALIZED", + "1,1,1", "global", "verification"); + const auto metrics = openDataset( + file.get(), std::string{kStepRoot} + "/global/verification_metrics"); + EXPECT_EQ( + readStringAttribute(metrics.get(), "metric_definition_ids"), + "free-residual-l2-over-max-free-force-l2,force-balance-l2-over-max-force-sum,moment-balance-l2-over-max-moment-sum"); + EXPECT_EQ( + readStringAttribute(metrics.get(), "acceptance_thresholds"), + "1e-10,1e-10,1e-10"); +} + +// MITC4-H5-003 +TEST(Hdf5ResultsWriter, WritesShellInventoryDespiteRequestsAndOmitsForbiddenPaths) { + TempDirectory directory{"shell-mandatory"}; + auto fixture = makeShellFixture(directory.path() / "shell.inp"); + const fesa::Diagnostic ignoredRequest{ + fesa::Severity::warning, + "ignored-output-request", + {fixture.domain->sourcePath(), 80U}, + "*ELEMENT OUTPUT", + "S", + "Output requests cannot filter mandatory shell results."}; + const auto output = directory.path() / "results.h5"; + + fesa::Hdf5ResultsWriter writer; + ASSERT_TRUE( + writer.write(output, *fixture.domain, *fixture.state, {ignoredRequest}) + .isOk()); + const auto file = openFile(output); + for (const char* suffix : { + "/element/shell/local_frame", + "/element/shell/generalized_strain", + "/element/shell/section_resultant", + "/element/shell/stress", + "/global/energy", + "/global/equilibrium", + "/global/verification_metrics"}) { + const std::string path = std::string{kStepRoot} + suffix; + EXPECT_GT(H5Lexists(file.get(), path.c_str(), H5P_DEFAULT), 0) << path; + } + for (const char* forbidden : { + "/steps/Step-1/frames/0/element/shell/drilling", + "/steps/Step-1/frames/0/element/shell/drilling_energy", + "/steps/Step-1/frames/0/element/shell/S33", + "/steps/Step-1/frames/0/element/shell/S13", + "/steps/Step-1/frames/0/element/shell/S23"}) { + EXPECT_EQ(H5Lexists(file.get(), forbidden, H5P_DEFAULT), 0) << forbidden; + } + EXPECT_EQ(datasetDimensions(file.get(), "/diagnostics"), + std::vector({1U})); +} + +// MITC4-H5-004 +TEST(Hdf5ResultsWriter, InvalidShellInventoryPreservesExistingFinal) { + TempDirectory directory{"shell-atomic"}; + auto fixture = makeShellFixture(directory.path() / "shell.inp"); + auto invalidState = fesa::AnalysisState::create( + *fixture.dofs, {"Step-1", 0U}); + const auto final = directory.path() / "results.h5"; + const std::vector sentinel = {'s', 'h', 'e', 'l', 'l'}; + writeBytes(final, sentinel); + + fesa::Hdf5ResultsWriter writer; + expectOutputFailure( + writer.write(final, *fixture.domain, invalidState, {}), + "invalid-result-rows"); + EXPECT_EQ(readBytes(final), sentinel); + EXPECT_EQ(entryCount(directory.path()), 1U); +}