From 5fe57cb145400a5e948078282bc7d5b744c9792f Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sun, 2 Aug 2026 01:14:40 +0900 Subject: [PATCH] =?UTF-8?q?feat(result-contract-completion):=20step=200=20?= =?UTF-8?q?=E2=80=94=20beam-element-end-recovery?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- include/fesa/elements/beam/beam3d2.hpp | 17 ++ phases/index.json | 2 +- src/fesa/assembly/parallel_assembler.cpp | 1 + src/fesa/assembly/serial_assembler.cpp | 1 + src/fesa/elements/beam/beam3d2.cpp | 125 +++++++++++++-- tests/CMakeLists.txt | 18 +++ tests/unit/elements/beam3d2_test.cpp | 188 +++++++++++++++++++++++ 7 files changed, 341 insertions(+), 11 deletions(-) diff --git a/include/fesa/elements/beam/beam3d2.hpp b/include/fesa/elements/beam/beam3d2.hpp index 265f8ce..268d2f5 100644 --- a/include/fesa/elements/beam/beam3d2.hpp +++ b/include/fesa/elements/beam/beam3d2.hpp @@ -2,22 +2,34 @@ #include #include +#include #include #include #include #include #include +#include #include namespace fesa { struct Beam3D2Input final { std::array coordinates; + std::array node_ids; IsotropicElastic material; BeamSection section; }; +struct BeamSectionResult final { + double xi; + NodeId end_node; + std::array section_strain; + std::array section_force; + double centroid_sigma_xx; + std::vector sigma_xx; +}; + struct Beam3D2Contribution final { Matrix12 local_stiffness; Matrix12 global_stiffness; @@ -32,4 +44,9 @@ struct BeamKernelResult final { [[nodiscard]] BeamKernelResult compute_beam3d2( const Beam3D2Input& input); +[[nodiscard]] std::vector recover_beam3d2( + const Beam3D2Input& input, + std::span element_displacement, + std::span> recovery_points); + } // namespace fesa diff --git a/phases/index.json b/phases/index.json index 2581dc9..200b69f 100644 --- a/phases/index.json +++ b/phases/index.json @@ -48,4 +48,4 @@ "status": "pending" } ] -} \ No newline at end of file +} diff --git a/src/fesa/assembly/parallel_assembler.cpp b/src/fesa/assembly/parallel_assembler.cpp index c7285e8..4f6ddd7 100644 --- a/src/fesa/assembly/parallel_assembler.cpp +++ b/src/fesa/assembly/parallel_assembler.cpp @@ -81,6 +81,7 @@ ElementEvaluation evaluate_element( domain.node(element.nodes[0]).position, domain.node(element.nodes[1]).position, }, + element.nodes, domain.material(element.material), domain.section(element.section), }); diff --git a/src/fesa/assembly/serial_assembler.cpp b/src/fesa/assembly/serial_assembler.cpp index 0904da2..35ac65f 100644 --- a/src/fesa/assembly/serial_assembler.cpp +++ b/src/fesa/assembly/serial_assembler.cpp @@ -128,6 +128,7 @@ std::vector collect_numeric_contributions( domain.node(element.nodes[0]).position, domain.node(element.nodes[1]).position, }, + element.nodes, domain.material(element.material), domain.section(element.section), }); diff --git a/src/fesa/elements/beam/beam3d2.cpp b/src/fesa/elements/beam/beam3d2.cpp index 5b98082..8b418c3 100644 --- a/src/fesa/elements/beam/beam3d2.cpp +++ b/src/fesa/elements/beam/beam3d2.cpp @@ -4,8 +4,10 @@ #include #include #include +#include #include #include +#include #include #include @@ -32,6 +34,20 @@ bool is_positive_finite(const double value) { return std::isfinite(value) && value > 0.0; } +std::array constitutive_values( + const Beam3D2Input& input) { + const double shear_modulus = + input.material.young / (2.0 * (1.0 + input.material.poisson)); + return { + input.material.young * input.section.area, + shear_modulus * input.section.shear_area_y, + shear_modulus * input.section.shear_area_z, + shear_modulus * input.section.torsion_j, + input.material.young * input.section.iy, + input.material.young * input.section.iz, + }; +} + std::optional validate_properties( const Beam3D2Input& input) { if (!std::isfinite(input.material.young) || @@ -149,6 +165,34 @@ Matrix12 transform_stiffness( return global; } +std::array transform_displacement( + const Matrix12& transformation, + const std::span global_displacement) { + std::array local_displacement{}; + for (std::size_t row = 0; row < transformation.size(); ++row) { + for (std::size_t column = 0; + column < transformation[row].size(); + ++column) { + local_displacement[row] += + transformation[row][column] * global_displacement[column]; + } + } + return local_displacement; +} + +std::array evaluate_strain( + const StrainMatrix& strain_matrix, + const std::array& local_displacement) { + std::array strain{}; + for (std::size_t component = 0; component < strain.size(); ++component) { + for (std::size_t dof = 0; dof < local_displacement.size(); ++dof) { + strain[component] += + strain_matrix[component][dof] * local_displacement[dof]; + } + } + return strain; +} + bool is_finite(const Matrix12& matrix) { for (const auto& row : matrix) { for (const double value : row) { @@ -189,16 +233,8 @@ BeamKernelResult compute_beam3d2(const Beam3D2Input& input) { "Beam kernel requires a representable positive Jacobian."); } - const double shear_modulus = - input.material.young / (2.0 * (1.0 + input.material.poisson)); - const std::array constitutive{ - input.material.young * input.section.area, - shear_modulus * input.section.shear_area_y, - shear_modulus * input.section.shear_area_z, - shear_modulus * input.section.torsion_j, - input.material.young * input.section.iy, - input.material.young * input.section.iz, - }; + const std::array constitutive = + constitutive_values(input); for (const double value : constitutive) { if (!is_positive_finite(value)) { return error_result( @@ -241,4 +277,73 @@ BeamKernelResult compute_beam3d2(const Beam3D2Input& input) { }; } +std::vector recover_beam3d2( + const Beam3D2Input& input, + const std::span element_displacement, + const std::span> recovery_points) { + const BeamKernelResult kernel = compute_beam3d2(input); + if (!kernel.contribution.has_value()) { + const std::string message = kernel.diagnostics.empty() + ? "Beam recovery requires a valid Beam3D2 input." + : kernel.diagnostics.front().message; + throw std::invalid_argument{message}; + } + + const Vec3 axis{ + input.coordinates[1].x - input.coordinates[0].x, + input.coordinates[1].y - input.coordinates[0].y, + input.coordinates[1].z - input.coordinates[0].z, + }; + const double jacobian = std::hypot(axis.x, axis.y, axis.z) / 2.0; + const Matrix12 transformation = + beam_transformation(kernel.contribution->frame); + const std::array local_displacement = + transform_displacement(transformation, element_displacement); + const std::array center_strain = evaluate_strain( + strain_matrix(0.0, jacobian), + local_displacement); + const std::array constitutive = + constitutive_values(input); + + std::vector results; + results.reserve(input.node_ids.size()); + for (std::size_t end = 0; end < input.node_ids.size(); ++end) { + const double xi = end == 0 ? -1.0 : 1.0; + std::array section_strain = evaluate_strain( + strain_matrix(xi, jacobian), + local_displacement); + section_strain[1] = center_strain[1]; + section_strain[2] = center_strain[2]; + + std::array section_force{}; + for (std::size_t component = 0; + component < section_force.size(); + ++component) { + section_force[component] = + constitutive[component] * section_strain[component]; + } + + std::vector sigma_xx; + sigma_xx.reserve(recovery_points.size()); + for (const auto& point : recovery_points) { + const double y = point[0]; + const double z = point[1]; + sigma_xx.push_back( + input.material.young * + (section_strain[0] + z * section_strain[4] - + y * section_strain[5])); + } + + results.push_back({ + xi, + input.node_ids[end], + section_strain, + section_force, + section_force[0] / input.section.area, + std::move(sigma_xx), + }); + } + return results; +} + } // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ca3f6b2..548c040 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -360,6 +360,24 @@ add_test( --gtest_filter=RigidBody.* ) +add_test( + NAME BeamRecovery + COMMAND "$" + --gtest_filter=BeamRecovery.* +) + +add_test( + NAME SectionForce + COMMAND "$" + --gtest_filter=SectionForce.* +) + +add_test( + NAME CentroidStress + COMMAND "$" + --gtest_filter=CentroidStress.* +) + add_executable(fesa_serial_assembly_tests unit/assembly/serial_assembler_test.cpp ) diff --git a/tests/unit/elements/beam3d2_test.cpp b/tests/unit/elements/beam3d2_test.cpp index 98b2dd5..a04ff46 100644 --- a/tests/unit/elements/beam3d2_test.cpp +++ b/tests/unit/elements/beam3d2_test.cpp @@ -22,6 +22,7 @@ constexpr double kShearAreaZ = 0.2; fesa::Beam3D2Input make_x_axis_input(const double length) { return { {{{0.0, 0.0, 0.0}, {length, 0.0, 0.0}}}, + {fesa::NodeId{1}, fesa::NodeId{2}}, {fesa::MaterialId{1}, "elastic", kYoung, kPoisson}, { fesa::SectionId{1}, @@ -369,4 +370,191 @@ TEST(Beam3D2, PreservesGlobalEnergyUnderRigidCoordinateRotation) { 4096.0 * std::numeric_limits::epsilon()); } +TEST(BeamRecovery, RecoversPureAxialStrainAndForceAtBothEnds) { + fesa::Beam3D2Input input = make_x_axis_input(2.0); + input.node_ids = {fesa::NodeId{101}, fesa::NodeId{202}}; + std::array displacement{}; + displacement[6] = 0.4; + + const auto results = fesa::recover_beam3d2(input, displacement, {}); + + ASSERT_EQ(results.size(), 2U); + EXPECT_DOUBLE_EQ(results[0].xi, -1.0); + EXPECT_EQ(results[0].end_node, fesa::NodeId{101}); + EXPECT_DOUBLE_EQ(results[1].xi, 1.0); + EXPECT_EQ(results[1].end_node, fesa::NodeId{202}); + for (const auto& result : results) { + expect_relative_near(result.section_strain[0], 0.2); + expect_relative_near(result.section_force[0], 16.8); + expect_relative_near(result.centroid_sigma_xx, 42.0); + for (std::size_t component = 1; component < 6; ++component) { + expect_relative_near(result.section_strain[component], 0.0); + expect_relative_near(result.section_force[component], 0.0); + } + EXPECT_TRUE(result.sigma_xx.empty()); + } +} + +TEST(BeamRecovery, RecoversPureTorsionAtBothEnds) { + fesa::Beam3D2Input input = make_x_axis_input(2.0); + input.node_ids = {fesa::NodeId{101}, fesa::NodeId{202}}; + std::array displacement{}; + displacement[9] = 0.6; + const double shear_modulus = kYoung / (2.0 * (1.0 + kPoisson)); + + const auto results = fesa::recover_beam3d2(input, displacement, {}); + + ASSERT_EQ(results.size(), 2U); + for (const auto& result : results) { + expect_relative_near(result.section_strain[3], 0.3); + expect_relative_near( + result.section_force[3], + shear_modulus * kTorsionJ * 0.3); + for (const std::size_t component : + std::array{0, 1, 2, 4, 5}) { + expect_relative_near(result.section_strain[component], 0.0); + expect_relative_near(result.section_force[component], 0.0); + } + } +} + +TEST(BeamRecovery, UsesReducedIntegrationPointForShearAtBothEnds) { + fesa::Beam3D2Input input = make_x_axis_input(2.0); + std::array displacement{}; + displacement[7] = 1.2; + displacement[5] = 0.2; + displacement[11] = 0.6; + displacement[8] = -0.4; + displacement[4] = 0.1; + displacement[10] = 0.5; + const double shear_modulus = kYoung / (2.0 * (1.0 + kPoisson)); + + const auto results = fesa::recover_beam3d2(input, displacement, {}); + + ASSERT_EQ(results.size(), 2U); + for (const auto& result : results) { + expect_relative_near(result.section_strain[1], 0.2); + expect_relative_near(result.section_strain[2], 0.1); + expect_relative_near( + result.section_force[1], + shear_modulus * kShearAreaY * 0.2); + expect_relative_near( + result.section_force[2], + shear_modulus * kShearAreaZ * 0.1); + } +} + +TEST(BeamRecovery, UsesBeamFrameForGlobalDisplacements) { + fesa::Beam3D2Input input = make_x_axis_input(2.0); + input.coordinates = {{{0.0, 0.0, 0.0}, {0.0, 2.0, 0.0}}}; + input.section.orientation = {0.0, 0.0, 1.0}; + input.node_ids = {fesa::NodeId{101}, fesa::NodeId{202}}; + std::array displacement{}; + displacement[7] = 0.4; + displacement[10] = 0.6; + const double shear_modulus = kYoung / (2.0 * (1.0 + kPoisson)); + + const auto results = fesa::recover_beam3d2(input, displacement, {}); + + ASSERT_EQ(results.size(), 2U); + for (const auto& result : results) { + expect_relative_near(result.section_strain[0], 0.2); + expect_relative_near(result.section_force[0], 16.8); + expect_relative_near(result.section_strain[3], 0.3); + expect_relative_near( + result.section_force[3], + shear_modulus * kTorsionJ * 0.3); + } +} + +TEST(SectionForce, RecoversPureBendingAboutLocalYAtBothEnds) { + fesa::Beam3D2Input input = make_x_axis_input(2.0); + input.node_ids = {fesa::NodeId{101}, fesa::NodeId{202}}; + std::array displacement{}; + displacement[4] = -0.4; + displacement[10] = 0.4; + + const auto results = fesa::recover_beam3d2(input, displacement, {}); + + ASSERT_EQ(results.size(), 2U); + for (const auto& result : results) { + expect_relative_near(result.section_strain[4], 0.4); + expect_relative_near(result.section_force[4], 2.52); + for (const std::size_t component : + std::array{0, 1, 2, 3, 5}) { + expect_relative_near(result.section_strain[component], 0.0); + expect_relative_near(result.section_force[component], 0.0); + } + } +} + +TEST(SectionForce, RecoversPureBendingAboutLocalZAtBothEnds) { + fesa::Beam3D2Input input = make_x_axis_input(2.0); + input.node_ids = {fesa::NodeId{101}, fesa::NodeId{202}}; + std::array displacement{}; + displacement[5] = 0.25; + displacement[11] = -0.25; + + const auto results = fesa::recover_beam3d2(input, displacement, {}); + + ASSERT_EQ(results.size(), 2U); + for (const auto& result : results) { + expect_relative_near(result.section_strain[5], -0.25); + expect_relative_near(result.section_force[5], -2.625); + for (const std::size_t component : + std::array{0, 1, 2, 3, 4}) { + expect_relative_near(result.section_strain[component], 0.0); + expect_relative_near(result.section_force[component], 0.0); + } + } +} + +TEST(SectionForce, PreservesBiaxialBendingSignsAtBothEnds) { + fesa::Beam3D2Input input = make_x_axis_input(2.0); + input.node_ids = {fesa::NodeId{101}, fesa::NodeId{202}}; + std::array displacement{}; + displacement[4] = -0.4; + displacement[10] = 0.4; + displacement[5] = 0.25; + displacement[11] = -0.25; + + const auto results = fesa::recover_beam3d2(input, displacement, {}); + + ASSERT_EQ(results.size(), 2U); + for (const auto& result : results) { + expect_relative_near(result.section_force[4], 2.52); + expect_relative_near(result.section_force[5], -2.625); + } +} + +TEST(CentroidStress, UsesAxialStressAndPreservesRecoveryPointOrder) { + fesa::Beam3D2Input input = make_x_axis_input(2.0); + input.node_ids = {fesa::NodeId{101}, fesa::NodeId{202}}; + std::array displacement{}; + displacement[6] = 0.2; + displacement[4] = -0.2; + displacement[10] = 0.2; + displacement[5] = 0.3; + displacement[11] = -0.3; + const std::array, 3> recovery_points{{ + {2.0, 3.0}, + {-1.0, 4.0}, + {5.0, -2.0}, + }}; + + const auto results = fesa::recover_beam3d2( + input, + displacement, + recovery_points); + + ASSERT_EQ(results.size(), 2U); + for (const auto& result : results) { + expect_relative_near(result.centroid_sigma_xx, 21.0); + ASSERT_EQ(result.sigma_xx.size(), recovery_points.size()); + expect_relative_near(result.sigma_xx[0], 273.0); + expect_relative_near(result.sigma_xx[1], 126.0); + expect_relative_near(result.sigma_xx[2], 252.0); + } +} + } // namespace