From 68253148459db1368bcfd3e842809a552e152255 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Fri, 9 Oct 2026 01:42:33 +0900 Subject: [PATCH] feat: bind Solid definitions to Hex8 runtime elements --- docs/hex8-structural-analysis/formulation.md | 12 +++ include/fesa/elements/element_factory.h | 4 + src/fesa/CMakeLists.txt | 2 + src/fesa/elements/element_factory.cpp | 69 +++++++++++++++ tests/CMakeLists.txt | 2 + tests/unit/elements/element_factory_test.cpp | 93 ++++++++++++++++++++ 6 files changed, 182 insertions(+) diff --git a/docs/hex8-structural-analysis/formulation.md b/docs/hex8-structural-analysis/formulation.md index 83ae6cd..0320733 100644 --- a/docs/hex8-structural-analysis/formulation.md +++ b/docs/hex8-structural-analysis/formulation.md @@ -154,10 +154,22 @@ Ma+f−F의 constrained projection이다. Initial u/v0에서 Mff a0=Ff(0+)−ff( Physical momenta는 interpolated current x/v를 사용해 \(p=\int\rho v\,dV_0\), \(L=\int x\times\rho v\,dV_0\)다.8점 consistent quadrature와 nodal mass 표현이 일치해야 한다. Angular-momentum rate는 \(\int x\times\rho a\,dV_0\)이며 v×v 항은0이다. +Consistent symmetric M에서 nodal 표현의 \(\sum v_a\times(Mv)_a=0\)도 성립한다. External/support moment는 같은 current 위치로 구한다. HHT 및 force trapezoid work는 frame별 exact energy conservation을 의미하지 않는다. Wext는 endpoint force·du trapezoid 누적, error=T+U−T0−U0−Wext다. Alpha0 linear unforced benchmark에서는 독립 보존을 확인한다. +Current moment의 검사에는 선형화의 적용 범위를 구분한다. Full physical residual을 +\(r=Ma+f_{int}-F\), 공개 constrained reaction을 \(Q=\Pi_c r\), free residual을 +\(e=\Pi_f r\)라 하면 \(B_x=\sum x_a\times(F_a+Q_a)-\dot L +=\sum x_a\times f_{int,a}-\sum x_a\times e_a\)다. 정적에서는 관성항을0으로 둔다. +TL은 유한회전 불변성으로 \(\sum x\times f_{int}=0\)지만 선형 K0는 +\(\sum X\times K_0u=0\)이며 \(\sum x\times K_0u=\sum u\times K_0u\)가 일반적으로 +O(u²)이다. 선형 current 수지 차이를 exact0으로 보정하거나 roundoff gate로 거부하지 않는다. +Current moment와 momentum은 그대로 계산하고 이 차이와 reference-arm identity를 내부 +evidence/로그로 설명한다. Newton 허용오차의 free residual 기여도 함께 남긴다. 시간차분 +\(\Delta L/h\)를 instantaneous \(\dot L\)와 exact equality로 판정하지 않는다. + ## 정식화 검증 경계 Analytic brick/affine energy, independent tensor polynomial integrals,6 rigid modes, scaling/ diff --git a/include/fesa/elements/element_factory.h b/include/fesa/elements/element_factory.h index b07b077..c33a308 100644 --- a/include/fesa/elements/element_factory.h +++ b/include/fesa/elements/element_factory.h @@ -29,6 +29,10 @@ class ElementFactory { /// @brief Creates one checked MITC4 runtime candidate. Result> CreateShell( const ElementDefinition& definition, const Domain& domain) const; + + /// @brief Creates one checked three-component C3D8 runtime candidate. + Result> CreateSolid( + const ElementDefinition& definition, const Domain& domain) const; }; } // namespace fesa diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index 67f5fcc..bdd099c 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -35,6 +35,8 @@ add_library( core/diagnostic.cpp core/status.cpp elements/element_factory.cpp + elements/hex8_element.cpp + elements/hex8_geometry.cpp elements/euler_beam_3d.cpp elements/finite_rotation.cpp elements/mitc4_shell.cpp diff --git a/src/fesa/elements/element_factory.cpp b/src/fesa/elements/element_factory.cpp index 52fec17..5fffbfa 100644 --- a/src/fesa/elements/element_factory.cpp +++ b/src/fesa/elements/element_factory.cpp @@ -9,6 +9,7 @@ #include #include "fesa/elements/euler_beam_3d.h" +#include "fesa/elements/hex8_element.h" #include "fesa/elements/mitc4_shell.h" #include "fesa/materials/material.h" #include "fesa/model/domain.h" @@ -225,6 +226,72 @@ Result> ElementFactory::CreateShell( return Result>::Success(std::move(result)); } +Result> ElementFactory::CreateSolid( + const ElementDefinition& definition, const Domain& domain) const { + const auto element_index = FindElementIndex(definition, domain); + const Hex8Definition* solid_definition = nullptr; + for (std::size_t index = 0U; index < domain.SolidElements().Size(); ++index) { + const auto& solid = domain.SolidElements()[index]; + if (&solid == &definition) { + solid_definition = &solid; + break; + } + } + if (!element_index.has_value() || solid_definition == nullptr) { + return FactoryFailure( + definition, domain, "invalid-element-definition", + "The C3D8 definition must be owned by the supplied Domain."); + } + if (!HasValidNodes(definition, domain, 8U)) { + return FactoryFailure(definition, domain, "invalid-element-topology", + "C3D8 requires eight valid Domain node indices."); + } + if (definition.PropertyIndex() >= domain.Properties().Size() || + definition.MaterialIndex() >= domain.Materials().Size()) { + return FactoryFailure( + definition, domain, "invalid-element-assignment", + "C3D8 property and material indices must resolve in the Domain."); + } + const auto& property = domain.Properties()[definition.PropertyIndex()]; + const auto& material = domain.Materials()[definition.MaterialIndex()]; + if (property.Kind() != ElementPropertyKind::kSolidSection || + material.Kind() != MaterialKind::kIsotropicLinearElastic) { + return FactoryFailure( + definition, domain, "incompatible-element-assignment", + "C3D8 requires a Solid section and isotropic elastic material."); + } + const auto& section = static_cast(property); + if (section.MaterialIndex() != definition.MaterialIndex()) { + return FactoryFailure( + definition, domain, "incompatible-element-assignment", + "C3D8 section and definition must select the same material."); + } + Hex8Element::Coordinates coordinates{}; + for (std::size_t node = 0U; node < coordinates.size(); ++node) { + coordinates[node] = domain.Nodes()[definition.NodeIndices()[node]].coordinates; + } + std::optional density; + for (const auto& record : domain.Densities()) { + if (record.material_index == definition.MaterialIndex()) { + density = record.rho; + break; + } + } + ElementDofLayout layout{definition.SourceId(), definition.NodeIndices(), + {DofComponent::kUx, DofComponent::kUy, + DofComponent::kUz}}; + auto candidate = Hex8Element::Create( + *element_index, std::move(layout), coordinates, + static_cast(material), density, + solid_definition->location); + if (!candidate.HasValue()) { + return Result>::Failure(candidate.GetStatus()); + } + std::unique_ptr result = + std::make_unique(std::move(candidate.Value())); + return Result>::Success(std::move(result)); +} + Result> ElementFactory::Create( const ElementDefinition& definition, const Domain& domain) const { switch (definition.Kind()) { @@ -232,6 +299,8 @@ Result> ElementFactory::Create( return CreateBeam(definition, domain); case ElementDefinitionKind::kMitc4Shell: return CreateShell(definition, domain); + case ElementDefinitionKind::kHex8: + return CreateSolid(definition, domain); } return FactoryFailure(definition, domain, "unsupported-element-definition", "The element definition kind is not supported."); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 25871db..0a81cf1 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -34,6 +34,7 @@ add_executable( unit/core/source_identity_test.cpp unit/core/status_test.cpp unit/elements/element_factory_test.cpp + unit/elements/hex8_element_test.cpp unit/elements/euler_beam_3d_test.cpp unit/elements/finite_rotation_test.cpp unit/elements/mitc4_shell_test.cpp @@ -50,6 +51,7 @@ add_executable( unit/math/vector_test.cpp unit/io/abaqus/domain_mapping_components_test.cpp unit/io/abaqus/domain_mapper_test.cpp + unit/io/abaqus/hex8_semantics_test.cpp unit/io/abaqus/eigen_domain_test.cpp unit/io/abaqus/hht_dynamic_domain_test.cpp unit/io/abaqus/thermal_domain_test.cpp diff --git a/tests/unit/elements/element_factory_test.cpp b/tests/unit/elements/element_factory_test.cpp index e41d97e..a69b960 100644 --- a/tests/unit/elements/element_factory_test.cpp +++ b/tests/unit/elements/element_factory_test.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -16,6 +17,7 @@ #include "fesa/elements/element.h" #include "fesa/elements/euler_beam_3d.h" +#include "fesa/elements/hex8_element.h" #include "fesa/elements/mitc4_shell.h" #include "fesa/fem/dof_manager.h" #include "fesa/model/domain.h" @@ -90,6 +92,38 @@ fesa::Domain MakeDomain(fesa::ModelDefinition definition) { return std::move(candidate.Value()); } +fesa::ModelDefinition SolidDefinition(bool with_density = true) { + fesa::ModelDefinition definition{}; + definition.source_path = "element-factory-solid.inp"; + definition.source_content_identity = "fnv1a64:solid"; + const fesa::Hex8Element::Coordinates coordinates{{ + {0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {1.0, 1.0, 0.0}, {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}, {1.0, 0.0, 1.0}, {1.0, 1.0, 1.0}, {0.0, 1.0, 1.0}}}; + for (std::size_t node = 0U; node < coordinates.size(); ++node) { + const auto label = static_cast(10U * (node + 1U)); + definition.nodes.push_back({{"Solid-1", label, std::to_string(label)}, + coordinates[node], + {definition.source_path, node + 1U}}); + } + definition.materials = {{"SolidMaterial", 1000.0, 0.3, + {definition.source_path, 10U}}}; + definition.solid_sections.emplace_back( + fesa::SourceEntityId{"Solid-1", 31, "31"}, "SolidSet", 0U, + fesa::SourceLocation{definition.source_path, 11U}); + definition.solid_elements.emplace_back( + fesa::SourceEntityId{"Solid-1", 99, "99"}, + std::array{0U, 1U, 2U, 3U, 4U, 5U, 6U, 7U}, + 0U, 0U, fesa::SourceLocation{definition.source_path, 12U}); + if (with_density) { + definition.densities.push_back({0U, 2.0, {definition.source_path, 13U}}); + } + definition.steps.emplace_back( + "Solid", std::vector{}, + std::vector{}, 1.0, 1.0, 1.0, 1.0, + fesa::SourceLocation{definition.source_path, 14U}); + return definition; +} + void ExpectSameMatrix(const fesa::Matrix& actual, const fesa::Matrix& expected) { ASSERT_EQ(actual.Rows(), expected.Rows()); @@ -344,3 +378,62 @@ TEST(ElementFactory, RejectsUnknownAndIncompatibleDefinitions) { EXPECT_EQ(missing_frame_result.GetStatus().Diagnostics()[0U].code, "missing-shell-frame"); } + +TEST(ElementFactory, SolidBindingPreservesThreeComponentsAndReferenceMass) { + const auto domain = MakeDomain(SolidDefinition()); + auto candidate = fesa::ElementFactory{}.Create(domain.Elements()[0U], domain); + ASSERT_TRUE(candidate.HasValue()); + const auto* solid = dynamic_cast(candidate.Value().get()); + ASSERT_NE(solid, nullptr); + EXPECT_EQ(solid->DofLayout().source_id.source_label, 99); + EXPECT_EQ(solid->DofLayout().components_per_node, + (std::vector{fesa::DofComponent::kUx, + fesa::DofComponent::kUy, + fesa::DofComponent::kUz})); + auto stiffness = solid->ComputeStiffness(); + ASSERT_TRUE(stiffness.HasValue()); + EXPECT_EQ(stiffness.Value().values.Rows(), 24U); + auto mass = solid->ComputeMass(); + ASSERT_TRUE(mass.HasValue()); + double total_mass = 0.0; + for (std::size_t a = 0U; a < 8U; ++a) { + for (std::size_t b = 0U; b < 8U; ++b) { + total_mass += mass.Value().values(3U * a, 3U * b); + } + } + EXPECT_NEAR(total_mass, 2.0, 1.0e-13); + EXPECT_NEAR(mass.Value().values(0U, 3U), 2.0 / 54.0, 1.0e-15); + auto response = solid->EvaluateLinear(fesa::Vector{24U}); + ASSERT_TRUE(response.HasValue()); + for (std::size_t point = 0U; point < 8U; ++point) { + EXPECT_EQ(response.Value().points[point].element, 0U); + EXPECT_EQ(response.Value().points[point].integration_point, point + 1U); + } +} + +TEST(ElementFactory, SolidStaticBindingKeepsAbsentDensityAndRejectsForeignOwner) { + const auto domain = MakeDomain(SolidDefinition(false)); + auto candidate = fesa::ElementFactory{}.Create(domain.Elements()[0U], domain); + ASSERT_TRUE(candidate.HasValue()); + EXPECT_TRUE(candidate.Value()->ComputeStiffness().HasValue()); + const auto* solid = dynamic_cast(candidate.Value().get()); + ASSERT_NE(solid, nullptr); + EXPECT_FALSE(solid->ComputeMass().HasValue()); + const auto foreign = SolidDefinition(false); + auto rejected = fesa::ElementFactory{}.Create(foreign.solid_elements[0U], domain); + ASSERT_FALSE(rejected.HasValue()); + EXPECT_EQ(rejected.GetStatus().Diagnostics()[0U].code, + "invalid-element-definition"); +} + +TEST(ElementFactory, SolidBindingRejectsReferenceInversionWithSourceContext) { + auto definition = SolidDefinition(); + for (auto& node : definition.nodes) { + node.coordinates[0U] = -node.coordinates[0U]; + } + const auto domain = MakeDomain(std::move(definition)); + auto rejected = fesa::ElementFactory{}.Create(domain.Elements()[0U], domain); + ASSERT_FALSE(rejected.HasValue()); + EXPECT_EQ(rejected.GetStatus().Diagnostics()[0U].entity_identity, "Solid-1:99"); + EXPECT_EQ(rejected.GetStatus().Diagnostics()[0U].location.line, 12U); +}