diff --git a/src/fesa/io/abaqus/domain_mapper.cpp b/src/fesa/io/abaqus/domain_mapper.cpp index af70de9..9357958 100644 --- a/src/fesa/io/abaqus/domain_mapper.cpp +++ b/src/fesa/io/abaqus/domain_mapper.cpp @@ -14,6 +14,7 @@ #include #include +#include "fesa/math/vector3.h" #include "fesa/model/shell_geometry.h" namespace fesa { @@ -1800,61 +1801,41 @@ class MappingContext { bool ValidateGeometry(const RawElement& raw, const Node& first, const Node& second, const GeneralBeamSection& section) { - const auto norm = [](const std::array& vector) { - return std::hypot(vector[0], vector[1], vector[2]); - }; - const auto maximum_absolute = [](const std::array& vector) { + const auto maximum_absolute = [](const Vector3& vector) { return std::max( {std::abs(vector[0]), std::abs(vector[1]), std::abs(vector[2])}); }; + const Vector3 first_position{first.coordinates}; + const Vector3 second_position{second.coordinates}; // Compare both approved inequalities after a common scaling. This // preserves the exact ratios while avoiding overflow in x*x and in // subtraction between large finite coordinates. const double global_coordinate_scale = - std::max({1.0, maximum_absolute(first.coordinates), - maximum_absolute(second.coordinates)}); - std::array first_scaled{}; - std::array second_scaled{}; - std::array delta_scaled{}; - for (std::size_t coordinate = 0U; coordinate < 3U; ++coordinate) { - first_scaled[coordinate] = - first.coordinates[coordinate] / global_coordinate_scale; - second_scaled[coordinate] = - second.coordinates[coordinate] / global_coordinate_scale; - delta_scaled[coordinate] = - second_scaled[coordinate] - first_scaled[coordinate]; - } - const double length_ratio = norm(delta_scaled); + std::max({1.0, maximum_absolute(first_position), + maximum_absolute(second_position)}); + const Vector3 first_scaled = first_position / global_coordinate_scale; + const Vector3 second_scaled = second_position / global_coordinate_scale; + const Vector3 delta_scaled = second_scaled - first_scaled; + const double length_ratio = delta_scaled.Norm(); const double coordinate_norm_ratio = - std::max({1.0 / global_coordinate_scale, norm(first_scaled), - norm(second_scaled)}); + std::max({1.0 / global_coordinate_scale, first_scaled.Norm(), + second_scaled.Norm()}); if (!(length_ratio > 1.0e-12 * coordinate_norm_ratio)) { return ModelFailure( "invalid-beam-length", raw.location, "ELEMENT", raw.label_text, "Beam length fails the approved scale-aware threshold."); } - std::array tangent{delta_scaled[0] / length_ratio, - delta_scaled[1] / length_ratio, - delta_scaled[2] / length_ratio}; + const Vector3 tangent = delta_scaled / length_ratio; - const double global_guide_scale = - std::max(1.0, maximum_absolute(section.first_axis)); - std::array guide_scaled{}; - for (std::size_t coordinate = 0U; coordinate < 3U; ++coordinate) { - guide_scaled[coordinate] = - section.first_axis[coordinate] / global_guide_scale; - } - const double projection = guide_scaled[0] * tangent[0] + - guide_scaled[1] * tangent[1] + - guide_scaled[2] * tangent[2]; - std::array perpendicular{ - guide_scaled[0] - projection * tangent[0], - guide_scaled[1] - projection * tangent[1], - guide_scaled[2] - projection * tangent[2]}; + const Vector3 guide{section.first_axis}; + const double global_guide_scale = std::max(1.0, maximum_absolute(guide)); + const Vector3 guide_scaled = guide / global_guide_scale; + const double projection = guide_scaled.Dot(tangent); + const Vector3 perpendicular = guide_scaled - projection * tangent; const double guide_norm_ratio = - std::max(1.0 / global_guide_scale, norm(guide_scaled)); - if (!(norm(perpendicular) > 1.0e-12 * guide_norm_ratio)) { + std::max(1.0 / global_guide_scale, guide_scaled.Norm()); + if (!(perpendicular.Norm() > 1.0e-12 * guide_norm_ratio)) { return ModelFailure( "invalid-beam-guide-vector", section.location, "BEAM GENERAL SECTION", raw.label_text, diff --git a/src/fesa/io/hdf5/hdf5_results_writer.cpp b/src/fesa/io/hdf5/hdf5_results_writer.cpp index c04794b..e912b2a 100644 --- a/src/fesa/io/hdf5/hdf5_results_writer.cpp +++ b/src/fesa/io/hdf5/hdf5_results_writer.cpp @@ -20,6 +20,7 @@ #include "fesa/analysis/analysis_model.h" #include "fesa/build_info.h" #include "fesa/fem/dof_manager.h" +#include "fesa/math/vector3.h" namespace fesa { namespace { @@ -163,11 +164,6 @@ Status OutputFailure(const std::string& code, const std::string& message) { {{Severity::kError, code, {}, "", "", message}}); } -bool IsFinite(const std::array& values) { - return std::all_of(values.begin(), values.end(), - [](const double value) { return std::isfinite(value); }); -} - template bool IsFinite(const std::array& values) { return std::all_of(values.begin(), values.end(), @@ -244,26 +240,21 @@ const char* ShellSourceTypeName(const ShellSourceElementType type) { bool IsOrthonormalRightHanded( const std::array, 3>& frame) { constexpr double kTolerance = 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) > kTolerance) { + const std::array axes = {Vector3{frame[0U]}, Vector3{frame[1U]}, + Vector3{frame[2U]}}; + for (const Vector3& axis : axes) { + if (!axis.IsFinite() || std::abs(axis.Dot(axis) - 1.0) > kTolerance) { return false; } } - if (std::abs(dot(frame[0U], frame[1U])) > kTolerance || - std::abs(dot(frame[0U], frame[2U])) > kTolerance || - std::abs(dot(frame[1U], frame[2U])) > kTolerance) { + if (std::abs(axes[0U].Dot(axes[1U])) > kTolerance || + std::abs(axes[0U].Dot(axes[2U])) > kTolerance || + std::abs(axes[1U].Dot(axes[2U])) > kTolerance) { 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) <= kTolerance; + const Vector3 cross = axes[0U].Cross(axes[1U]); + const double handedness = cross.Dot(axes[2U]); + return handedness > 0.0 && std::abs(handedness - 1.0) <= kTolerance; } bool ComputeLocalAxes(const Domain& domain, @@ -275,30 +266,25 @@ bool ComputeLocalAxes(const Domain& domain, } const auto& first = domain.Nodes()[element.node_indices[0U]].coordinates; const auto& second = domain.Nodes()[element.node_indices[1U]].coordinates; - const auto& guide = domain.Sections()[element.section_index].first_axis; - const std::array delta = { - second[0U] - first[0U], second[1U] - first[1U], second[2U] - first[2U]}; - const double length = std::hypot(delta[0U], delta[1U], delta[2U]); - if (!IsFinite(first) || !IsFinite(second) || !IsFinite(guide) || - !IsFinite(delta) || !std::isfinite(length) || !(length > 0.0)) { + const Vector3 first_position{first}; + const Vector3 second_position{second}; + const Vector3 guide{domain.Sections()[element.section_index].first_axis}; + const Vector3 delta = second_position - first_position; + const double length = delta.Norm(); + if (!first_position.IsFinite() || !second_position.IsFinite() || + !guide.IsFinite() || !delta.IsFinite() || !std::isfinite(length) || + !(length > 0.0)) { return false; } - const std::array x = {delta[0U] / length, delta[1U] / length, - delta[2U] / length}; - const double projection = - guide[0U] * x[0U] + guide[1U] * x[1U] + guide[2U] * x[2U]; - const std::array y_trial = {guide[0U] - projection * x[0U], - guide[1U] - projection * x[1U], - guide[2U] - projection * x[2U]}; - const double y_norm = std::hypot(y_trial[0U], y_trial[1U], y_trial[2U]); - if (!IsFinite(y_trial) || !std::isfinite(y_norm) || !(y_norm > 0.0)) { + const Vector3 x = delta / length; + const double projection = guide.Dot(x); + const Vector3 y_trial = guide - projection * x; + const double y_norm = y_trial.Norm(); + if (!y_trial.IsFinite() || !std::isfinite(y_norm) || !(y_norm > 0.0)) { return false; } - const std::array y = {y_trial[0U] / y_norm, y_trial[1U] / y_norm, - y_trial[2U] / y_norm}; - const std::array z = {x[1U] * y[2U] - x[2U] * y[1U], - x[2U] * y[0U] - x[0U] * y[2U], - x[0U] * y[1U] - x[1U] * y[0U]}; + const Vector3 y = y_trial / y_norm; + const Vector3 z = x.Cross(y); axes = {x[0U], x[1U], x[2U], y[0U], y[1U], y[2U], z[0U], z[1U], z[2U]}; return IsFinite(axes); } @@ -417,7 +403,7 @@ Status ValidateShellWriterInput(const Domain& domain, ++position) { if (row.stress[position].position != positions[position] || row.stress[position].zeta != kZeta[position] || - !IsFinite(row.stress[position].components)) { + !Vector3{row.stress[position].components}.IsFinite()) { return OutputFailure("invalid-result-rows", "Shell stress rows must preserve BOTTOM, " "MIDDLE, TOP identity."); @@ -426,7 +412,7 @@ Status ValidateShellWriterInput(const Domain& domain, } if (!std::isfinite(state.PhysicalStrainEnergy()) || !IsFinite(state.Equilibrium()) || - !IsFinite(state.VerificationMetrics())) { + !Vector3{state.VerificationMetrics()}.IsFinite()) { return OutputFailure( "invalid-result-rows", "Shell energy, equilibrium, and verification metrics must be finite."); @@ -493,7 +479,7 @@ Status ValidateWriterInput(const std::filesystem::path& output_path, node.source_id.source_label_text.empty() || !IsValidUtf8(node.source_id.instance_name) || !IsValidUtf8(node.source_id.source_label_text) || - !IsFinite(node.coordinates)) { + !Vector3{node.coordinates}.IsFinite()) { return OutputFailure( "invalid-result-identity", "Every node requires finite coordinates and UTF-8 source identity."); diff --git a/src/fesa/results/result_recovery.cpp b/src/fesa/results/result_recovery.cpp index 69cb79f..d3401f3 100644 --- a/src/fesa/results/result_recovery.cpp +++ b/src/fesa/results/result_recovery.cpp @@ -16,6 +16,7 @@ #include "fesa/elements/euler_beam_3d.h" #include "fesa/elements/mitc4_shell.h" +#include "fesa/math/vector3.h" namespace fesa { namespace { @@ -389,27 +390,12 @@ Result> ResolveLoadTarget(const Domain& domain, "A station-eligibility load target must resolve to a node or node set."); } -std::array Cross(const std::array& left, - const std::array& right) { - return {left[1U] * right[2U] - left[2U] * right[1U], - left[2U] * right[0U] - left[0U] * right[2U], - left[0U] * right[1U] - left[1U] * right[0U]}; -} - -double Dot(const std::array& left, - const std::array& right) { - return left[0U] * right[0U] + left[1U] * right[1U] + left[2U] * right[2U]; -} - -double Norm(const std::array& value) { - return std::hypot(value[0U], value[1U], value[2U]); -} - bool AccumulateVectorAndScale(std::array& total, double& scale, const std::array& contribution) { - const double magnitude = Norm(contribution); + const Vector3 contribution_vector{contribution}; + const double magnitude = contribution_vector.Norm(); const double accumulated_scale = scale + magnitude; - if (!IsFinite(contribution) || !std::isfinite(magnitude) || + if (!contribution_vector.IsFinite() || !std::isfinite(magnitude) || !std::isfinite(accumulated_scale)) { return false; } @@ -427,7 +413,7 @@ bool AccumulateVectorAndScale(std::array& total, double& scale, double NormalizedBalance(const std::array& balance, const double scale) { - const double balance_norm = Norm(balance); + const double balance_norm = Vector3{balance}.Norm(); if (!std::isfinite(balance_norm) || !std::isfinite(scale)) { return (std::numeric_limits::infinity)(); } @@ -474,10 +460,11 @@ Status PopulateShellGlobalEvidence(const Domain& domain, const DofManager& dofs, } } - const auto applied_force_moment = - Cross(domain.Nodes()[node].coordinates, nodal_applied_force); - const auto reaction_force_moment = - Cross(domain.Nodes()[node].coordinates, nodal_reaction_force); + const Vector3 position{domain.Nodes()[node].coordinates}; + const Vector3 applied_force_moment = + position.Cross(Vector3{nodal_applied_force}); + const Vector3 reaction_force_moment = + position.Cross(Vector3{nodal_reaction_force}); for (std::size_t component = 0U; component < 3U; ++component) { nodal_applied_moment[component] += applied_force_moment[component]; nodal_reaction_moment[component] += reaction_force_moment[component]; @@ -498,24 +485,24 @@ Status PopulateShellGlobalEvidence(const Domain& domain, const DofManager& dofs, } } - std::array force_balance{}; - std::array moment_balance{}; + const Vector3 force_balance = + Vector3{applied_force} + Vector3{reaction_force}; + const Vector3 moment_balance = + Vector3{applied_moment} + Vector3{reaction_moment}; for (std::size_t component = 0U; component < 3U; ++component) { - force_balance[component] = - applied_force[component] + reaction_force[component]; - moment_balance[component] = - applied_moment[component] + reaction_moment[component]; candidate.equilibrium[component] = force_balance[component]; candidate.equilibrium[3U + component] = moment_balance[component]; } - const double force_metric = NormalizedBalance( - force_balance, (std::max)(applied_force_scale, reaction_force_scale)); + const double force_metric = + NormalizedBalance(force_balance.Components(), + (std::max)(applied_force_scale, reaction_force_scale)); const double moment_metric = NormalizedBalance( - moment_balance, (std::max)(applied_moment_scale, reaction_moment_scale)); + moment_balance.Components(), + (std::max)(applied_moment_scale, reaction_moment_scale)); candidate.verification_metrics = {normalized_residual, force_metric, moment_metric}; if (!IsFinite(candidate.equilibrium) || - !IsFinite(candidate.verification_metrics)) { + !Vector3{candidate.verification_metrics}.IsFinite()) { return RecoveryFailure("nonfinite-recovery-value", {domain.SourcePath(), 0U}, "global-equilibrium", "Global equilibrium values and their physical " @@ -535,32 +522,25 @@ std::optional LocalAxes(const Domain& domain, const EulerBeam3DDefinition& element) { const auto& first = domain.Nodes()[element.node_indices[0U]].coordinates; const auto& second = domain.Nodes()[element.node_indices[1U]].coordinates; - const std::array delta = { - second[0U] - first[0U], second[1U] - first[1U], second[2U] - first[2U]}; - const double length = Norm(delta); + const Vector3 delta = Vector3{second} - Vector3{first}; + const double length = delta.Norm(); if (!std::isfinite(length) || !(length > 0.0)) { return std::nullopt; } - const std::array ex = {delta[0U] / length, delta[1U] / length, - delta[2U] / length}; - const auto& guide = domain.Sections()[element.section_index].first_axis; - const double projection = Dot(guide, ex); - const std::array ey_trial = {guide[0U] - projection * ex[0U], - guide[1U] - projection * ex[1U], - guide[2U] - projection * ex[2U]}; - const double ey_norm = Norm(ey_trial); + const Vector3 ex = delta / length; + const Vector3 guide{domain.Sections()[element.section_index].first_axis}; + const double projection = guide.Dot(ex); + const Vector3 ey_trial = guide - projection * ex; + const double ey_norm = ey_trial.Norm(); if (!std::isfinite(ey_norm) || !(ey_norm > 0.0)) { return std::nullopt; } - const std::array ey = { - ey_trial[0U] / ey_norm, ey_trial[1U] / ey_norm, ey_trial[2U] / ey_norm}; - const std::array ez = Cross(ex, ey); - const AxisSet axes = {ex, ey, ez}; + const Vector3 ey = ey_trial / ey_norm; + const Vector3 ez = ex.Cross(ey); + const AxisSet axes = {ex.Components(), ey.Components(), ez.Components()}; for (const auto& axis : axes) { - for (const double component : axis) { - if (!std::isfinite(component)) { - return std::nullopt; - } + if (!Vector3{axis}.IsFinite()) { + return std::nullopt; } } return axes; diff --git a/tests/unit/io/abaqus/domain_mapper_test.cpp b/tests/unit/io/abaqus/domain_mapper_test.cpp index 546042b..d8863b4 100644 --- a/tests/unit/io/abaqus/domain_mapper_test.cpp +++ b/tests/unit/io/abaqus/domain_mapper_test.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -13,6 +14,7 @@ #include #include "fesa/io/abaqus/input_reader.h" +#include "fesa/math/vector3.h" namespace { @@ -239,6 +241,7 @@ RootAssembly, 6, -12.5 } // namespace +// C-DUP-002 TEST(InpDomainMapping, MapsEverySupportedKeywordAndLegacyDeck) { auto result = MapText("supported-inventory", SupportedInventoryDeck(true)); ASSERT_TRUE(result.HasValue()); @@ -425,6 +428,17 @@ OnlySecond, 2, 5. "2, 1., 0., 0.", "2, 1e308, 1e297, 0."), "0., 1., 0.", "1e308, 0., 0.")); ASSERT_TRUE(large_finite.HasValue()); + const auto& large_domain = large_finite.Value(); + ASSERT_EQ(large_domain.Nodes().size(), 2U); + const fesa::Vector3 large_delta = + fesa::Vector3{large_domain.Nodes()[1U].coordinates} - + fesa::Vector3{large_domain.Nodes()[0U].coordinates}; + EXPECT_TRUE(large_delta.IsFinite()); + EXPECT_TRUE(std::isfinite(large_delta.Norm())); + EXPECT_DOUBLE_EQ(large_delta.Y(), 1.0e297); + const fesa::Vector3 large_guide{large_domain.Sections()[0U].first_axis}; + EXPECT_TRUE(large_guide.IsFinite()); + EXPECT_TRUE(std::isfinite(large_guide.Norm())); } TEST(InpDomainMapping, KeepsNodeAndElementSetNamesInSeparateNamespaces) { diff --git a/tests/unit/io/hdf5/hdf5_results_writer_test.cpp b/tests/unit/io/hdf5/hdf5_results_writer_test.cpp index ab1e5d5..2cef031 100644 --- a/tests/unit/io/hdf5/hdf5_results_writer_test.cpp +++ b/tests/unit/io/hdf5/hdf5_results_writer_test.cpp @@ -24,6 +24,7 @@ #include "fesa/analysis/analysis_state.h" #include "fesa/build_info.h" #include "fesa/fem/dof_manager.h" +#include "fesa/math/vector3.h" #include "fesa/model/domain.h" namespace { @@ -1114,7 +1115,7 @@ TEST(Hdf5ResultsWriter, SuccessfullyReplacesExistingFinal) { EXPECT_EQ(ReadUint64Attribute(metadata.Get(), "schema_version"), 0U); } -// MITC4-H5-001 +// MITC4-H5-001, C-DUP-002 TEST(Hdf5ResultsWriter, WritesExactShellMetadataAndModelIdentity) { TempDirectory directory{"shell-model"}; const auto source = directory.Path() / "shell.inp"; @@ -1150,9 +1151,15 @@ TEST(Hdf5ResultsWriter, WritesExactShellMetadataAndModelIdentity) { EXPECT_EQ(DatasetDimensions(file.Get(), "/model/shell/nodal_director"), std::vector({4U, 3U})); + constexpr fesa::Vector3 director{0.0, 0.0, 1.0}; + std::vector expected_directors; + for (std::size_t node = 0U; node < 4U; ++node) { + expected_directors.insert(expected_directors.end(), + director.Components().begin(), + director.Components().end()); + } 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})); + expected_directors); EXPECT_EQ(DatasetDimensions(file.Get(), "/model/shell/nodal_frame"), std::vector({4U, 3U, 3U})); ExpectCompoundMemberNames(file.Get(), "/model/shell/materials", @@ -1305,3 +1312,31 @@ TEST(Hdf5ResultsWriter, InvalidShellInventoryPreservesExistingFinal) { EXPECT_EQ(ReadBytes(final), sentinel); EXPECT_EQ(EntryCount(directory.Path()), 1U); } + +// C-DUP-002 +TEST(Hdf5ResultsWriter, + Vector3ShellSerializationRejectsNonfiniteWithoutReplacingFinal) { + TempDirectory directory{"shell-vector3-atomic"}; + auto fixture = MakeShellFixture(directory.Path() / "shell.inp"); + const auto final = directory.Path() / "results.h5"; + fesa::Hdf5ResultsWriter writer; + + ASSERT_TRUE(writer.Write(final, *fixture.domain, *fixture.state, {}).IsOk()); + const std::vector valid_bytes = ReadBytes(final); + const auto file = OpenFile(final); + const auto directors = + ReadDoubleDataset(file.Get(), "/model/shell/nodal_director"); + ASSERT_EQ(directors.size(), 12U); + for (std::size_t offset = 0U; offset < directors.size(); offset += 3U) { + const fesa::Vector3 director{directors[offset], directors[offset + 1U], + directors[offset + 2U]}; + EXPECT_TRUE(director.IsFinite()); + EXPECT_EQ(director, (fesa::Vector3{0.0, 0.0, 1.0})); + } + + fixture.state->Displacement()[0U] = std::numeric_limits::quiet_NaN(); + ExpectOutputFailure(writer.Write(final, *fixture.domain, *fixture.state, {}), + "invalid-result-state"); + EXPECT_EQ(ReadBytes(final), valid_bytes); + EXPECT_EQ(EntryCount(directory.Path()), 1U); +} diff --git a/tests/unit/results/result_recovery_test.cpp b/tests/unit/results/result_recovery_test.cpp index 487eb60..0ba1260 100644 --- a/tests/unit/results/result_recovery_test.cpp +++ b/tests/unit/results/result_recovery_test.cpp @@ -20,6 +20,7 @@ #include "fesa/assembly/parallel_for.h" #include "fesa/assembly/sparse_assembler.h" #include "fesa/fem/dof_manager.h" +#include "fesa/math/vector3.h" #include "fesa/model/domain.h" namespace { @@ -719,14 +720,14 @@ TEST(ResultRecovery, KeepsFullResidualAndComputesGlobalShellEquilibrium) { } } -// MITC4-REC-004 +// MITC4-REC-004, C-DUP-002 TEST(ResultRecovery, UsesGlobalOriginForShellMomentBalance) { auto centered_definition = MakeShellDefinition(); auto translated_definition = centered_definition; - constexpr std::array translation{7.0, 11.0, 0.0}; + constexpr fesa::Vector3 translation{7.0, 11.0, 0.0}; for (auto& node : translated_definition.nodes) { - for (std::size_t component = 0U; component < translation.size(); - ++component) { + for (std::size_t component = 0U; + component < translation.Components().size(); ++component) { node.coordinates[component] += translation[component]; } } @@ -753,13 +754,8 @@ TEST(ResultRecovery, UsesGlobalOriginForShellMomentBalance) { EXPECT_NEAR(translated.Equilibrium()[component], centered_force[component], 1.0e-12); } - const std::array translated_moment_delta{ - translation[1U] * centered_force[2U] - - translation[2U] * centered_force[1U], - translation[2U] * centered_force[0U] - - translation[0U] * centered_force[2U], - translation[0U] * centered_force[1U] - - translation[1U] * centered_force[0U]}; + const fesa::Vector3 translated_moment_delta = + translation.Cross(fesa::Vector3{centered_force}); for (std::size_t component = 0U; component < 3U; ++component) { EXPECT_NEAR(translated.Equilibrium()[3U + component] - centered.Equilibrium()[3U + component],