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

test: add mitc4 patch benchmark coverage

Browse Source
This commit is contained in:
NINI
2026-05-04 23:05:01 +09:00
parent 5465ea61d8
commit ebdb2519dc
5 changed files with 252 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/test_main.cpp
+207
View File
@@ -160,6 +160,76 @@ fesa::Real quadraticEnergy(const fesa::DenseMatrix& stiffness, const std::vector
return energy;
}
std::array<fesa::Vec3, 4> flatUnitSquareCoordinates() {
return {fesa::Vec3{0.0, 0.0, 0.0},
fesa::Vec3{1.0, 0.0, 0.0},
fesa::Vec3{1.0, 1.0, 0.0},
fesa::Vec3{0.0, 1.0, 0.0}};
}
using MITC4DofField = std::function<fesa::Vec3(const fesa::Vec3&)>;
fesa::MITC4ElementDofVector elementDofsFromFields(const std::array<fesa::Vec3, 4>& coords,
MITC4DofField translation,
MITC4DofField rotation) {
fesa::MITC4ElementDofVector values{};
for (std::size_t node = 0; node < coords.size(); ++node) {
const auto displacement = translation(coords[node]);
const auto rotational_dof = rotation(coords[node]);
values[node * 6 + 0] = displacement.x;
values[node * 6 + 1] = displacement.y;
values[node * 6 + 2] = displacement.z;
values[node * 6 + 3] = rotational_dof.x;
values[node * 6 + 4] = rotational_dof.y;
values[node * 6 + 5] = rotational_dof.z;
}
return values;
}
std::vector<fesa::Real> toVector(const fesa::MITC4ElementDofVector& values) {
return {values.begin(), values.end()};
}
fesa::MITC4StrainVector localStrainAt(const fesa::MITC4Geometry& geometry,
const fesa::MITC4ElementDofVector& values,
fesa::Real xi,
fesa::Real eta,
fesa::Real zeta) {
const auto rows = fesa::mitc4TiedCovariantStrainRows(geometry, xi, eta, zeta);
FESA_CHECK(rows.ok());
const auto convected_strain = fesa::evaluateMITC4StrainRows(rows, values);
FESA_CHECK(convected_strain.ok());
const auto basis = fesa::computeMITC4IntegrationBasis(geometry, xi, eta, zeta);
FESA_CHECK(basis.ok());
const auto transform = fesa::mitc4CovariantToLocalStrainTransform(basis);
FESA_CHECK(transform.ok());
return fesa::multiplyMITC4MaterialMatrix(transform.matrix, convected_strain.values);
}
fesa::Real singleElementCantileverTipUz(fesa::Real thickness) {
fesa::Domain domain;
domain.nodes[1] = {1, {0.0, 0.0, 0.0}};
domain.nodes[2] = {2, {4.0, 0.0, 0.0}};
domain.nodes[3] = {3, {4.0, 1.0, 0.0}};
domain.nodes[4] = {4, {0.0, 1.0, 0.0}};
domain.elements[1] = {1, fesa::ElementType::MITC4, {1, 2, 3, 4}, "SHELLS"};
domain.node_sets["clamped"] = {"CLAMPED", {1, 4}};
domain.element_sets["shells"] = {"SHELLS", {1}};
domain.materials["shell_steel"] = {"SHELL_STEEL", 210000.0, 0.30};
domain.shell_sections.push_back({"SHELLS", "SHELL_STEEL", thickness});
domain.boundary_conditions.push_back({"CLAMPED", 1, 6, 0.0});
domain.loads.push_back({"2", 3, -0.5});
domain.loads.push_back({"3", 3, -0.5});
fesa::LinearStaticAnalysis analysis;
const auto result = analysis.run(domain);
FESA_CHECK(result.ok());
fesa::DofManager dofs(domain);
const auto node2_uz = result.state.u_full[static_cast<std::size_t>(dofs.fullIndex(2, fesa::Dof::UZ))];
const auto node3_uz = result.state.u_full[static_cast<std::size_t>(dofs.fullIndex(3, fesa::Dof::UZ))];
return 0.5 * (node2_uz + node3_uz);
}
fesa::Domain singleElementValidationDomain() {
fesa::Domain domain;
domain.nodes[1] = {1, {0, 0, 0}};
@@ -1316,6 +1386,143 @@ FESA_TEST(mitc4_rigid_translation_energy_is_zero_and_drilling_energy_is_document
FESA_CHECK_NEAR(quadraticEnergy(result.global, drilling), 4.0 * result.drilling_stiffness, 1.0e-8);
}
FESA_TEST(mitc4_constant_membrane_patch_strain_is_uniform) {
const auto coords = flatUnitSquareCoordinates();
const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
FESA_CHECK(geometry.ok());
const fesa::Real eps_x = 0.010;
const fesa::Real eps_y = -0.004;
const fesa::Real gamma_xy = 0.006;
const auto dofs = elementDofsFromFields(
coords,
[&](const fesa::Vec3& p) {
return fesa::Vec3{eps_x * p.x + 0.5 * gamma_xy * p.y,
eps_y * p.y + 0.5 * gamma_xy * p.x,
0.0};
},
[](const fesa::Vec3&) { return fesa::Vec3{}; });
for (const auto& sample : fesa::mitc4GaussQuadrature2x2x2()) {
const auto strain = localStrainAt(geometry, dofs, sample.xi, sample.eta, sample.zeta);
FESA_CHECK_NEAR(strain[0], eps_x, 1.0e-12);
FESA_CHECK_NEAR(strain[1], eps_y, 1.0e-12);
FESA_CHECK_NEAR(strain[2], 0.0, 1.0e-12);
FESA_CHECK_NEAR(strain[3], 0.0, 1.0e-12);
FESA_CHECK_NEAR(strain[4], 0.0, 1.0e-12);
FESA_CHECK_NEAR(strain[5], gamma_xy, 1.0e-12);
}
const auto stiffness = fesa::mitc4ElementStiffness(coords, 210000.0, 0.30, 0.2);
FESA_CHECK(stiffness.ok());
const auto energy = quadraticEnergy(stiffness.local_without_drilling, toVector(dofs));
FESA_CHECK(energy > 0.0);
}
FESA_TEST(mitc4_pure_bending_patch_is_membrane_free_and_thickness_antisymmetric) {
const auto coords = flatUnitSquareCoordinates();
const fesa::Real thickness = 0.2;
const auto geometry = fesa::buildMITC4Geometry(coords, thickness);
FESA_CHECK(geometry.ok());
const fesa::Real curvature_x = 0.025;
const auto dofs = elementDofsFromFields(
coords,
[](const fesa::Vec3&) { return fesa::Vec3{}; },
[&](const fesa::Vec3& p) { return fesa::Vec3{0.0, curvature_x * p.x, 0.0}; });
const fesa::Real xi = -0.30;
const fesa::Real eta = 0.45;
const fesa::Real zeta = 1.0 / std::sqrt(3.0);
const auto top = localStrainAt(geometry, dofs, xi, eta, zeta);
const auto mid = localStrainAt(geometry, dofs, xi, eta, 0.0);
const auto bottom = localStrainAt(geometry, dofs, xi, eta, -zeta);
FESA_CHECK(std::fabs(top[0]) > 1.0e-6);
FESA_CHECK_NEAR(top[0] + bottom[0], 0.0, 1.0e-12);
FESA_CHECK_NEAR(mid[0], 0.0, 1.0e-12);
FESA_CHECK_NEAR(top[1] + bottom[1], 0.0, 1.0e-12);
FESA_CHECK_NEAR(mid[1], 0.0, 1.0e-12);
const auto stiffness = fesa::mitc4ElementStiffness(coords, 210000.0, 0.30, thickness);
FESA_CHECK(stiffness.ok());
const auto energy = quadraticEnergy(stiffness.local_without_drilling, toVector(dofs));
FESA_CHECK(energy > 0.0);
}
FESA_TEST(mitc4_pure_transverse_shear_patch_is_constant) {
const auto coords = flatUnitSquareCoordinates();
const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
FESA_CHECK(geometry.ok());
const fesa::Real gamma_xz = 0.012;
const auto dofs = elementDofsFromFields(
coords,
[&](const fesa::Vec3& p) { return fesa::Vec3{0.0, 0.0, gamma_xz * p.x}; },
[](const fesa::Vec3&) { return fesa::Vec3{}; });
for (const auto& sample : fesa::mitc4GaussQuadrature2x2x2()) {
const auto strain = localStrainAt(geometry, dofs, sample.xi, sample.eta, sample.zeta);
FESA_CHECK_NEAR(strain[0], 0.0, 1.0e-12);
FESA_CHECK_NEAR(strain[1], 0.0, 1.0e-12);
FESA_CHECK_NEAR(strain[2], 0.0, 1.0e-12);
FESA_CHECK_NEAR(strain[3], 0.0, 1.0e-12);
FESA_CHECK_NEAR(strain[4], gamma_xz, 1.0e-12);
FESA_CHECK_NEAR(strain[5], 0.0, 1.0e-12);
}
}
FESA_TEST(mitc4_pure_twist_patch_reproduces_bilinear_transverse_shear) {
const auto coords = flatUnitSquareCoordinates();
const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
FESA_CHECK(geometry.ok());
const fesa::Real twist = 0.018;
const auto dofs = elementDofsFromFields(
coords,
[&](const fesa::Vec3& p) { return fesa::Vec3{0.0, 0.0, twist * p.x * p.y}; },
[](const fesa::Vec3&) { return fesa::Vec3{}; });
for (const auto& sample : fesa::mitc4GaussQuadrature2x2x2()) {
const auto x = 0.5 * (sample.xi + 1.0);
const auto y = 0.5 * (sample.eta + 1.0);
const auto strain = localStrainAt(geometry, dofs, sample.xi, sample.eta, sample.zeta);
FESA_CHECK_NEAR(strain[0], 0.0, 1.0e-12);
FESA_CHECK_NEAR(strain[1], 0.0, 1.0e-12);
FESA_CHECK_NEAR(strain[2], 0.0, 1.0e-12);
FESA_CHECK_NEAR(strain[3], twist * x, 1.0e-12);
FESA_CHECK_NEAR(strain[4], twist * y, 1.0e-12);
FESA_CHECK_NEAR(strain[5], 0.0, 1.0e-12);
}
}
FESA_TEST(mitc4_membrane_patch_energy_is_not_hidden_by_drilling_scale) {
const auto coords = flatUnitSquareCoordinates();
const auto dofs = elementDofsFromFields(
coords,
[](const fesa::Vec3& p) { return fesa::Vec3{0.003 * p.x, -0.002 * p.y, 0.0}; },
[](const fesa::Vec3&) { return fesa::Vec3{}; });
fesa::ElementStiffnessOptions no_drilling;
no_drilling.drilling_stiffness_scale = 0.0;
fesa::ElementStiffnessOptions strong_drilling;
strong_drilling.drilling_stiffness_scale = 1.0e-1;
const auto physical = fesa::mitc4ElementStiffness(coords, 210000.0, 0.30, 0.2, no_drilling);
const auto stabilized = fesa::mitc4ElementStiffness(coords, 210000.0, 0.30, 0.2, strong_drilling);
FESA_CHECK(physical.ok());
FESA_CHECK(stabilized.ok());
const auto physical_energy = quadraticEnergy(physical.local_with_drilling, toVector(dofs));
const auto stabilized_energy = quadraticEnergy(stabilized.local_with_drilling, toVector(dofs));
FESA_CHECK(physical_energy > 0.0);
FESA_CHECK_NEAR(stabilized_energy, physical_energy, physical_energy * 1.0e-12);
}
FESA_TEST(mitc4_single_element_cantilever_is_thickness_sensitive) {
const auto thick_tip = singleElementCantileverTipUz(0.20);
const auto thin_tip = singleElementCantileverTipUz(0.10);
FESA_CHECK(thick_tip < 0.0);
FESA_CHECK(thin_tip < thick_tip);
const auto response_ratio = std::fabs(thin_tip / thick_tip);
FESA_CHECK(response_ratio > 2.0);
}
FESA_TEST(mitc4_internal_force_is_stiffness_times_displacement_for_linear_phase1) {
const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2);