baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat: bind Solid definitions to Hex8 runtime elements

Browse Source
This commit is contained in:
KOKO\Mimi
2026-10-09 01:42:33 +09:00
parent 744b34804d
commit 6825314845
6 changed files with 182 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/hex8-structural-analysis/formulation.md
+12
View File
@@ -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/
include/fesa/elements/element_factory.h
+4
View File
@@ -29,6 +29,10 @@ class ElementFactory {
/// @brief Creates one checked MITC4 runtime candidate.
Result<std::unique_ptr<Element>> CreateShell(
const ElementDefinition& definition, const Domain& domain) const;
/// @brief Creates one checked three-component C3D8 runtime candidate.
Result<std::unique_ptr<Element>> CreateSolid(
const ElementDefinition& definition, const Domain& domain) const;
};
} // namespace fesa
src/fesa/CMakeLists.txt
+2
View File
@@ -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
src/fesa/elements/element_factory.cpp
+69
View File
@@ -9,6 +9,7 @@
#include <vector>
#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<std::unique_ptr<Element>> ElementFactory::CreateShell(
return Result<std::unique_ptr<Element>>::Success(std::move(result));
}
Result<std::unique_ptr<Element>> 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<const SolidSection&>(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<double> 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<const IsotropicLinearElasticMaterial&>(material), density,
solid_definition->location);
if (!candidate.HasValue()) {
return Result<std::unique_ptr<Element>>::Failure(candidate.GetStatus());
}
std::unique_ptr<Element> result =
std::make_unique<Hex8Element>(std::move(candidate.Value()));
return Result<std::unique_ptr<Element>>::Success(std::move(result));
}
Result<std::unique_ptr<Element>> ElementFactory::Create(
const ElementDefinition& definition, const Domain& domain) const {
switch (definition.Kind()) {
@@ -232,6 +299,8 @@ Result<std::unique_ptr<Element>> 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.");
tests/CMakeLists.txt
+2
View File
@@ -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
tests/unit/elements/element_factory_test.cpp
+93
View File
@@ -8,6 +8,7 @@
#include <functional>
#include <memory>
#include <stdexcept>
#include <string>
#include <string_view>
#include <type_traits>
#include <utility>
@@ -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<std::int64_t>(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<fesa::EntityIndex, 8>{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<fesa::PrescribedDisplacementDefinition>{},
std::vector<fesa::NodalLoad>{}, 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<const fesa::Hex8Element*>(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>{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<const fesa::Hex8Element*>(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);
}