feat(linear-static-mitc4-shell): step 11 - shell-hdf5-output

This commit is contained in:
KOKO\Mimi
2026-08-12 21:42:14 +09:00
parent c060bb763b
commit dea6093cf2
2 changed files with 1517 additions and 99 deletions
@@ -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<fesa::EntityIndex>(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<fesa::Domain>(
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<fesa::DofManager>(
std::move(dofsResult.value()));
auto state = std::make_unique<fesa::AnalysisState>(
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<double>(index + 1U);
state->externalForce()[index] = 10.0 + static_cast<double>(index);
state->internalForce()[index] = 20.0 + static_cast<double>(index);
state->residual()[index] = 30.0 + static_cast<double>(index);
state->reaction()[index] = 40.0 + static_cast<double>(index);
}
const double gauss = 1.0 / std::sqrt(3.0);
const std::array<std::array<double, 2>, 4> coordinates{{
{-gauss, -gauss},
{gauss, -gauss},
{gauss, gauss},
{-gauss, gauss}}};
const std::array<fesa::ShellMidsurfaceLocation, 4> 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<double>(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<double> readDoubleDataset(
return values;
}
std::vector<std::uint8_t> 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<std::size_t>(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<std::uint8_t> 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<std::string>& 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<int>(expectedNames.size()));
for (std::size_t index = 0U; index < expectedNames.size(); ++index) {
char* rawName = H5Tget_member_name(
type.get(), static_cast<unsigned>(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<hsize_t>({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<hsize_t>({4U, 3U}));
EXPECT_EQ(
readDoubleDataset(file.get(), "/model/shell/nodal_director"),
std::vector<double>({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<hsize_t>({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<hsize_t>({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<double>({-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<hsize_t>({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<char> 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);
}