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

feat: add mitc4 covariant strain tying

Browse Source
This commit is contained in:
NINI
2026-05-04 22:02:25 +09:00
parent fa492f994d
commit d550e13e77
5 changed files with 369 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/test_main.cpp
+149
View File
@@ -1,5 +1,6 @@
#include "fesa/fesa.hpp"
#include <array>
#include <cmath>
#include <functional>
#include <iostream>
@@ -140,6 +141,12 @@ void checkRightHandedOrthonormal(const fesa::Vec3& e1, const fesa::Vec3& e2, con
FESA_CHECK_NEAR(fesa::dot(fesa::cross(e1, e2), e3), 1.0, 1.0e-14);
}
std::array<fesa::Real, 24> zeroElementDofs() {
std::array<fesa::Real, 24> values{};
values.fill(0.0);
return values;
}
fesa::Domain singleElementValidationDomain() {
fesa::Domain domain;
domain.nodes[1] = {1, {0, 0, 0}};
@@ -921,6 +928,148 @@ FESA_TEST(mitc4_geometry_reports_singular_geometry_and_thickness) {
FESA_CHECK(fesa::containsDiagnostic(corner_basis.diagnostics, "FESA-MITC4-SINGULAR-JACOBIAN"));
}
FESA_TEST(mitc4_rotation_transform_and_displacement_interpolation) {
const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
FESA_CHECK(geometry.ok());
const auto flat_rotation = fesa::mitc4LocalRotations(geometry.nodal_frames[0], {1.0, 2.0, 3.0});
FESA_CHECK_NEAR(flat_rotation.alpha, 1.0, 1.0e-14);
FESA_CHECK_NEAR(flat_rotation.beta, 2.0, 1.0e-14);
FESA_CHECK_NEAR(flat_rotation.gamma, 3.0, 1.0e-14);
checkVecNear(fesa::mitc4DirectorIncrement(geometry.nodal_frames[0], {0.0, 0.0, 5.0}), {0, 0, 0}, 1.0e-14);
checkVecNear(fesa::mitc4DirectorIncrement(geometry.nodal_frames[0], {0.0, 2.0, 5.0}), {2, 0, 0}, 1.0e-14);
std::array<fesa::Real, 24> dofs = zeroElementDofs();
for (std::size_t node = 0; node < 4; ++node) {
dofs[6 * node + 0] = 1.0;
dofs[6 * node + 1] = 0.5;
dofs[6 * node + 4] = 2.0;
dofs[6 * node + 5] = 99.0;
}
const auto mid = fesa::mitc4DisplacementDerivatives(geometry, dofs, 0.0, 0.0, 0.0);
FESA_CHECK(mid.ok());
checkVecNear(mid.displacement, {1.0, 0.5, 0.0}, 1.0e-14);
const auto top = fesa::mitc4DisplacementDerivatives(geometry, dofs, 0.0, 0.0, 1.0);
FESA_CHECK(top.ok());
checkVecNear(top.displacement, {1.2, 0.5, 0.0}, 1.0e-14);
const std::array<fesa::Vec3, 4> fallback_coords = {{{0, 0, 0}, {1, 0, 0}, {1, 0, -1}, {0, 0, -1}}};
const auto fallback_geometry = fesa::buildMITC4Geometry(fallback_coords, 0.1);
FESA_CHECK(fallback_geometry.ok());
const auto fallback_rotation = fesa::mitc4LocalRotations(fallback_geometry.nodal_frames[0], {2.0, 3.0, 5.0});
FESA_CHECK_NEAR(fallback_rotation.alpha, -2.0, 1.0e-14);
FESA_CHECK_NEAR(fallback_rotation.beta, 5.0, 1.0e-14);
FESA_CHECK_NEAR(fallback_rotation.gamma, 3.0, 1.0e-14);
}
FESA_TEST(mitc4_direct_covariant_strain_rows_match_finite_difference) {
const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
FESA_CHECK(geometry.ok());
constexpr fesa::Real xi = 0.2;
constexpr fesa::Real eta = -0.3;
constexpr fesa::Real zeta = 0.4;
const auto rows = fesa::mitc4DirectCovariantStrainRows(geometry, xi, eta, zeta);
FESA_CHECK(rows.ok());
FESA_CHECK(fesa::mitc4StrainComponentLabels()[0] == "eps11");
FESA_CHECK(fesa::mitc4StrainComponentLabels()[3] == "gamma23");
const std::array<std::size_t, 7> checked_dofs = {0, 1, 2, 7, 9, 10, 11};
const fesa::Real h = 1.0e-6;
for (std::size_t dof : checked_dofs) {
auto plus = zeroElementDofs();
auto minus = zeroElementDofs();
plus[dof] = h;
minus[dof] = -h;
const auto plus_strain = fesa::mitc4DirectCovariantStrain(geometry, plus, xi, eta, zeta);
const auto minus_strain = fesa::mitc4DirectCovariantStrain(geometry, minus, xi, eta, zeta);
FESA_CHECK(plus_strain.ok());
FESA_CHECK(minus_strain.ok());
for (std::size_t component = 0; component < 6; ++component) {
const fesa::Real finite_difference = (plus_strain.values[component] - minus_strain.values[component]) / (2.0 * h);
FESA_CHECK_NEAR(rows.rows[component][dof], finite_difference, 1.0e-8);
}
}
const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13);
const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23);
for (std::size_t node = 0; node < 4; ++node) {
const std::size_t drilling_dof = 6 * node + 5;
for (std::size_t component = 0; component < 6; ++component) {
FESA_CHECK_NEAR(rows.rows[component][drilling_dof], 0.0, 1.0e-14);
}
}
FESA_CHECK(gamma13 == 4);
FESA_CHECK(gamma23 == 3);
}
FESA_TEST(mitc4_mitc_tying_rows_use_fesa_ac_bd_signs) {
const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
FESA_CHECK(geometry.ok());
const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23);
const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13);
const auto direct_a = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, -1.0, 0.0);
const auto direct_b = fesa::mitc4DirectCovariantStrainRows(geometry, -1.0, 0.0, 0.0);
const auto direct_c = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, 1.0, 0.0);
const auto direct_d = fesa::mitc4DirectCovariantStrainRows(geometry, 1.0, 0.0, 0.0);
const auto mitc_a = fesa::mitc4TiedCovariantStrainRows(geometry, 0.0, -1.0, 0.0);
const auto mitc_b = fesa::mitc4TiedCovariantStrainRows(geometry, -1.0, 0.0, 0.0);
const auto mitc_c = fesa::mitc4TiedCovariantStrainRows(geometry, 0.0, 1.0, 0.0);
const auto mitc_d = fesa::mitc4TiedCovariantStrainRows(geometry, 1.0, 0.0, 0.0);
FESA_CHECK(direct_a.ok());
FESA_CHECK(direct_b.ok());
FESA_CHECK(direct_c.ok());
FESA_CHECK(direct_d.ok());
FESA_CHECK(mitc_a.ok());
FESA_CHECK(mitc_b.ok());
FESA_CHECK(mitc_c.ok());
FESA_CHECK(mitc_d.ok());
for (std::size_t dof = 0; dof < 24; ++dof) {
FESA_CHECK_NEAR(mitc_a.rows[gamma13][dof], direct_a.rows[gamma13][dof], 1.0e-14);
FESA_CHECK_NEAR(mitc_c.rows[gamma13][dof], direct_c.rows[gamma13][dof], 1.0e-14);
FESA_CHECK_NEAR(mitc_b.rows[gamma23][dof], direct_b.rows[gamma23][dof], 1.0e-14);
FESA_CHECK_NEAR(mitc_d.rows[gamma23][dof], direct_d.rows[gamma23][dof], 1.0e-14);
}
}
FESA_TEST(mitc4_mitc_gauss_shear_rows_are_interpolated_from_tying_rows) {
const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
FESA_CHECK(geometry.ok());
constexpr fesa::Real xi = 0.5;
constexpr fesa::Real eta = 0.25;
constexpr fesa::Real zeta = 0.0;
const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23);
const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13);
const auto direct_at_point = fesa::mitc4DirectCovariantStrainRows(geometry, xi, eta, zeta);
const auto direct_a = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, -1.0, zeta);
const auto direct_b = fesa::mitc4DirectCovariantStrainRows(geometry, -1.0, 0.0, zeta);
const auto direct_c = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, 1.0, zeta);
const auto direct_d = fesa::mitc4DirectCovariantStrainRows(geometry, 1.0, 0.0, zeta);
const auto tied = fesa::mitc4TiedCovariantStrainRows(geometry, xi, eta, zeta);
FESA_CHECK(direct_at_point.ok());
FESA_CHECK(tied.ok());
const std::size_t node2_ry = 6 + 4;
const fesa::Real expected_gamma13 =
0.5 * (1.0 - eta) * direct_a.rows[gamma13][node2_ry] + 0.5 * (1.0 + eta) * direct_c.rows[gamma13][node2_ry];
FESA_CHECK_NEAR(tied.rows[gamma13][node2_ry], expected_gamma13, 1.0e-14);
FESA_CHECK(std::fabs(tied.rows[gamma13][node2_ry] - direct_at_point.rows[gamma13][node2_ry]) > 1.0e-4);
const std::size_t node4_rx = 3 * 6 + 3;
const fesa::Real expected_gamma23 =
0.5 * (1.0 - xi) * direct_b.rows[gamma23][node4_rx] + 0.5 * (1.0 + xi) * direct_d.rows[gamma23][node4_rx];
FESA_CHECK_NEAR(tied.rows[gamma23][node4_rx], expected_gamma23, 1.0e-14);
FESA_CHECK(std::fabs(tied.rows[gamma23][node4_rx] - direct_at_point.rows[gamma23][node4_rx]) > 1.0e-4);
for (std::size_t component : {std::size_t{0}, std::size_t{1}, std::size_t{2}, std::size_t{5}}) {
for (std::size_t dof = 0; dof < 24; ++dof) {
FESA_CHECK_NEAR(tied.rows[component][dof], direct_at_point.rows[component][dof], 1.0e-14);
}
}
}
FESA_TEST(mitc4_shape_functions_and_stiffness_baseline) {
auto shape = fesa::shapeFunctions(0.25, -0.5);
const fesa::Real sum = shape.n[0] + shape.n[1] + shape.n[2] + shape.n[3];