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

feat: integrate transactional Solid static procedures

Browse Source
This commit is contained in:
KOKO\Mimi
2026-10-09 02:37:30 +09:00
parent 0a78d516eb a34466f7c3
commit 9ce2b26cd9
6 changed files with 1727 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/unit/analysis/solid_static_procedure_test.cpp
+826
View File
@@ -0,0 +1,826 @@
#include "fesa/analysis/solid_static_procedure.h"
#include <gtest/gtest.h>
#include <cmath>
#include <filesystem>
#include <fstream>
#include <iostream>
#include <limits>
#include <sstream>
#include <string>
#include "fesa/assembly/parallel_for.h"
#include "fesa/io/abaqus/domain_mapper.h"
#include "fesa/io/abaqus/input_reader.h"
#include "fesa/math/sparse_matrix.h"
#include "fesa/solvers/linear/mkl_pardiso_solver.h"
#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h"
namespace {
const double kX[8][3] = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0},
{0, 0, 1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}};
std::string Deck(const std::string& boundary, const std::string& load = "",
bool nonlinear = false, double h = 0.25, double period = 1.0) {
std::ostringstream out;
out << "*Part,name=P\n*Node\n";
for (int a = 0; a < 8; ++a)
out << a + 1 << ',' << kX[a][0] << ',' << kX[a][1] << ',' << kX[a][2]
<< '\n';
out << "*Element,type=C3D8,elset=E\n99,1,2,3,4,5,6,7,8\n"
"*Solid Section,elset=E,material=M\n,\n*End Part\n"
"*Assembly,name=A\n*Instance,name=I,part=P\n*End Instance\n"
"*End Assembly\n*Material,name=M\n*Elastic\n120,.25\n"
"*Step,name=Solid,inc=100";
if (nonlinear) out << ",nlgeom=YES";
out << "\n*Static,direct\n" << h << ',' << period << '\n';
if (!boundary.empty()) out << "*Boundary\n" << boundary;
if (!load.empty()) out << "*Cload\n" << load;
out << "*End Step\n";
return out.str();
}
fesa::Result<fesa::Domain> Map(const std::string& text) {
const auto path = std::filesystem::temp_directory_path() /
"fesa-solid-static-procedure.inp";
{
std::ofstream out{path};
out << text;
}
auto syntax = fesa::AbaqusInputReader{}.Read(path);
std::filesystem::remove(path);
if (!syntax.HasValue())
return fesa::Result<fesa::Domain>::Failure(syntax.GetStatus());
return fesa::AbaqusDomainMapper{}.Map(syntax.Value());
}
std::string Message(const fesa::Status& status) {
return status.Diagnostics().empty() ? "" : status.Diagnostics()[0].message;
}
fesa::Result<fesa::SolidStaticExecutionResult> ExecuteSolid(
const fesa::Domain& domain) {
fesa::SerialParallelFor serial;
fesa::MklPardisoSolver linear;
fesa::MklPardisoTangentSolver tangent;
return fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
.Execute(domain);
}
// Symmetry supports permit the exact affine uniaxial solution
// u=(e*x,-nu*e*y,-nu*e*z).
std::string Supports(double left = 0.0) {
std::ostringstream out;
for (int a = 0; a < 8; ++a) {
if (kX[a][0] == 0) out << a + 1 << ",1,1," << left << '\n';
if (kX[a][1] == 0) out << a + 1 << ",2,2,0\n";
if (kX[a][2] == 0) out << a + 1 << ",3,3,0\n";
}
return out.str();
}
std::string Tension(double force) {
std::ostringstream out;
for (int a = 0; a < 8; ++a)
if (kX[a][0] == 1) out << a + 1 << ",1," << force / 4.0 << '\n';
return out.str();
}
std::string AffineBoundary(double ex, double shear = 0.0, double angle = 0.0) {
std::ostringstream out;
out.precision(17);
for (int a = 0; a < 8; ++a) {
const double x = (1 + ex) * kX[a][0] + shear * kX[a][1];
const double y = kX[a][1];
const double u[3] = {std::cos(angle) * x - std::sin(angle) * y - kX[a][0],
std::sin(angle) * x + std::cos(angle) * y - kX[a][1],
0};
for (int i = 0; i < 3; ++i)
out << a + 1 << ',' << i + 1 << ',' << i + 1 << ',' << u[i] << '\n';
}
return out.str();
}
class FailedTangent final : public fesa::TangentSolver {
public:
enum class Mode { kBackend, kIteration, kNonfinite, kGeometry };
explicit FailedTangent(Mode mode, bool after_accepted = false)
: mode_(mode), after_accepted_(after_accepted) {}
fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override {
// First increment is a genuine accepted solve, then inject a backend fault.
++calls_;
if (mode_ == Mode::kBackend && calls_ > 3U)
return fesa::Status::Failure(fesa::FailureCategory::kSolver,
{{fesa::Severity::kError,
"injected-backend",
{},
"STATIC",
{},
"Backend failed."}});
return solver_.FactorizeFresh(matrix);
}
fesa::Status SolveCorrection(const fesa::Vector& rhs,
fesa::Vector& du) const override {
if (mode_ == Mode::kBackend || (after_accepted_ && calls_ <= 2U))
return solver_.SolveCorrection(rhs, du);
for (std::size_t i = 0; i < du.Size(); ++i)
du[i] = mode_ == Mode::kNonfinite
? std::numeric_limits<double>::quiet_NaN()
: mode_ == Mode::kGeometry ? -100.0
: 0.0;
return fesa::Status::Ok();
}
private:
Mode mode_;
bool after_accepted_;
unsigned calls_{0U};
fesa::MklPardisoTangentSolver solver_;
};
class ObservedIndefiniteTangent final : public fesa::TangentSolver {
public:
fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override {
rayleigh_sum = 0.0;
for (double value : matrix.Values()) rayleigh_sum += value;
return solver_.FactorizeFresh(matrix);
}
fesa::Status SolveCorrection(const fesa::Vector& rhs,
fesa::Vector& du) const override {
return solver_.SolveCorrection(rhs, du);
}
double rayleigh_sum{0.0};
private:
fesa::MklPardisoTangentSolver solver_;
};
} // namespace
TEST(SolidStaticProcedure, LinearForceNonzeroSupportAndCompleteRecovery) {
auto domain = Map(Deck(Supports(0.02), Tension(1.2)));
ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 1U);
const auto& frame = result.Value().state.AcceptedFrames()[0];
EXPECT_EQ(frame.frame_index, 0U);
EXPECT_DOUBLE_EQ(frame.load_factor, 1.0);
ASSERT_EQ(frame.displacement.Size(), 24U);
ASSERT_EQ(frame.solid_results.size(), 8U);
for (std::size_t a = 0; a < 8; ++a) {
EXPECT_NEAR(frame.displacement[3 * a], 0.02 + 0.01 * kX[a][0], 1e-13);
EXPECT_NEAR(frame.displacement[3 * a + 1], -0.0025 * kX[a][1], 1e-13);
EXPECT_NEAR(frame.displacement[3 * a + 2], -0.0025 * kX[a][2], 1e-13);
if (kX[a][0] == 1)
EXPECT_DOUBLE_EQ(frame.reaction[3 * a], 0.0);
else
EXPECT_NEAR(frame.reaction[3 * a], -0.3, 1e-13);
}
EXPECT_NEAR(frame.physical_strain_energy, 0.006, 1e-13);
for (const auto& point : frame.solid_results) {
EXPECT_NEAR(point.cauchy_stress[0], 1.2, 1e-12);
EXPECT_NEAR(point.cauchy_stress[1], 0.0, 1e-12);
EXPECT_NEAR(point.strain[0], 0.01, 1e-13);
}
for (double value : frame.global_equilibrium) EXPECT_NEAR(value, 0.0, 1e-12);
}
TEST(SolidStaticProcedure,
AllConstrainedShearHasTensorStrainAndLinearMomentEvidence) {
auto domain = Map(Deck(AffineBoundary(0.02, 0.03)));
ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
const auto& frame = result.Value().state.CommittedState();
for (const auto& point : frame.solid_results) {
EXPECT_NEAR(point.strain[3], 0.015, 1e-13);
EXPECT_NEAR(point.cauchy_stress[3], 1.44, 1e-12);
}
EXPECT_NEAR(frame.physical_strain_energy, 0.0504, 1e-13);
EXPECT_GT(std::abs(frame.balance.current_internal_moment[2]), 1e-5);
EXPECT_NEAR(frame.balance.current_internal_moment[2],
frame.balance.linearization_moment[2], 1e-12);
EXPECT_NEAR(frame.balance.reference_internal_moment[2], 0.0, 1e-12);
}
TEST(SolidStaticProcedure,
NonlinearDirectScheduleScalesForceAndAbsoluteTargets) {
auto domain = Map(Deck(Supports(0.02), Tension(1.2), true, 0.5, 2.0));
ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 4U);
for (std::size_t n = 0; n < 4; ++n) {
const auto& frame = result.Value().state.AcceptedFrames()[n];
EXPECT_EQ(frame.frame_index, n);
EXPECT_EQ(frame.increment_index, n + 1U);
EXPECT_DOUBLE_EQ(frame.load_factor, 0.25 * (n + 1));
EXPECT_DOUBLE_EQ(frame.displacement[0], 0.02 * frame.load_factor);
EXPECT_DOUBLE_EQ(frame.external_force[3], 0.3 * frame.load_factor);
EXPECT_GT(frame.iteration_count, 1U);
EXPECT_EQ(frame.is_factor_one, n == 3);
for (double value : frame.balance.current_internal_moment)
EXPECT_NEAR(value, 0.0, 1e-12);
}
ASSERT_FALSE(result.Value().state.IterationRecords().empty());
EXPECT_TRUE(result.Value().state.IterationRecords().back().accepted);
}
TEST(SolidStaticProcedure, PurePrescribedNonlinearAndSmallLoadLinearLimit) {
auto pure = Map(Deck(
Supports() + "2,1,1,.01\n3,1,1,.01\n6,1,1,.01\n7,1,1,.01\n", "", true));
ASSERT_TRUE(pure.HasValue());
auto result = ExecuteSolid(pure.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
for (const auto& record : result.Value().state.IterationRecords())
if (record.convergence_evaluated) {
EXPECT_EQ(record.force_residual.branch,
fesa::NewtonConvergenceBranch::kZeroScale);
EXPECT_DOUBLE_EQ(record.force_residual.denominator_scale, 0.0);
}
auto linear = Map(Deck(Supports(), Tension(1.2e-5)));
auto nonlinear = Map(Deck(Supports(), Tension(1.2e-5), true));
ASSERT_TRUE(linear.HasValue());
ASSERT_TRUE(nonlinear.HasValue());
auto l = ExecuteSolid(linear.Value());
auto n = ExecuteSolid(nonlinear.Value());
ASSERT_TRUE(l.HasValue());
ASSERT_TRUE(n.HasValue()) << Message(n.GetStatus());
for (std::size_t i = 0; i < 24; ++i)
EXPECT_NEAR(l.Value().state.CommittedState().displacement[i],
n.Value().state.CommittedState().displacement[i], 2e-13);
}
TEST(SolidStaticProcedure, StrainedThenRotatedTLCauchyAndEnergyOracle) {
const double angle = 1.2;
auto domain = Map(Deck(AffineBoundary(0.03, 0, angle), "", true));
ASSERT_TRUE(domain.HasValue());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
const auto& frame = result.Value().state.CommittedState();
const double green = (1.03 * 1.03 - 1) / 2;
const double s0 = 144 * green * 1.03;
const double s1 = 48 * green / 1.03;
for (const auto& point : frame.solid_results) {
EXPECT_NEAR(point.strain[0], green, 1e-13);
EXPECT_NEAR(
point.cauchy_stress[0],
std::pow(std::cos(angle), 2) * s0 + std::pow(std::sin(angle), 2) * s1,
1e-11);
EXPECT_NEAR(point.cauchy_stress[3],
std::sin(angle) * std::cos(angle) * (s0 - s1), 1e-11);
}
EXPECT_NEAR(frame.physical_strain_energy, 72 * green * green, 1e-12);
for (double value : frame.balance.current_internal_moment)
EXPECT_NEAR(value, 0.0, 1e-12);
}
TEST(SolidStaticProcedure, FailedBackendKeepsActualAcceptedHistory) {
auto domain = Map(Deck(Supports(), Tension(1.2), true));
ASSERT_TRUE(domain.HasValue());
fesa::SerialParallelFor serial;
fesa::MklPardisoSolver linear;
FailedTangent tangent{FailedTangent::Mode::kBackend};
fesa::SolidStaticState state;
auto status = fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
.ExecuteInto(domain.Value(), state);
ASSERT_FALSE(status.IsOk());
ASSERT_EQ(state.AcceptedFrames().size(), 1U);
EXPECT_DOUBLE_EQ(state.CommittedState().load_factor, 0.25);
EXPECT_EQ(state.CommittedState().displacement.Size(), 24U);
EXPECT_EQ(state.CommittedState().solid_results.size(), 8U);
EXPECT_EQ(state.IterationRecords().back().diagnostic_code,
"injected-backend");
auto bad = state.CommittedState();
bad.frame_index = 1;
bad.increment_index = 2;
bad.load_factor = 0.5;
bad.displacement[0] = std::numeric_limits<double>::quiet_NaN();
EXPECT_FALSE(state.CommitAcceptedFrame(bad).IsOk());
bad = state.CommittedState();
bad.frame_index = 1;
bad.increment_index = 2;
bad.load_factor = .5;
bad.solid_results.pop_back();
EXPECT_FALSE(state.CommitAcceptedFrame(bad).IsOk());
EXPECT_EQ(state.AcceptedFrames().size(), 1U);
EXPECT_DOUBLE_EQ(state.CommittedState().load_factor, 0.25);
}
TEST(SolidStaticProcedure,
GeometryIterationAndNonfiniteFailuresPreserveBaseline) {
auto domain = Map(Deck(Supports(), Tension(1.2), true));
ASSERT_TRUE(domain.HasValue());
for (auto mode :
{FailedTangent::Mode::kGeometry, FailedTangent::Mode::kIteration,
FailedTangent::Mode::kNonfinite}) {
fesa::SerialParallelFor serial;
fesa::MklPardisoSolver linear;
FailedTangent tangent{mode};
fesa::SolidStaticState state;
auto status = fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
.ExecuteInto(domain.Value(), state);
EXPECT_FALSE(status.IsOk());
EXPECT_TRUE(state.AcceptedFrames().empty());
EXPECT_DOUBLE_EQ(state.CommittedState().load_factor, 0.0);
ASSERT_EQ(state.CommittedState().displacement.Size(), 24U);
for (std::size_t i = 0; i < 24; ++i)
EXPECT_DOUBLE_EQ(state.CommittedState().displacement[i], 0.0);
ASSERT_FALSE(state.IterationRecords().empty());
EXPECT_FALSE(state.IterationRecords().back().accepted);
EXPECT_FALSE(state.IterationRecords().back().diagnostic_code.empty());
if (mode == FailedTangent::Mode::kIteration)
EXPECT_EQ(state.IterationRecords().back().iteration, 100U);
}
}
TEST(SolidStaticProcedure, UnconstrainedSolidFailsWithoutRegularization) {
auto domain = Map(Deck("", Tension(1.2)));
ASSERT_TRUE(domain.HasValue());
auto result = ExecuteSolid(domain.Value());
EXPECT_FALSE(result.HasValue());
}
TEST(SolidStaticProcedure,
EightElementAffinePatchHasSolvedInteriorAndStableRows) {
std::ostringstream out;
out << "*Part,name=P\n*Node\n";
const auto id = [](int x, int y, int z) { return 1 + x + 3 * y + 9 * z; };
for (int z = 0; z < 3; ++z)
for (int y = 0; y < 3; ++y)
for (int x = 0; x < 3; ++x)
out << id(x, y, z) << ',' << x << ',' << y << ',' << z << '\n';
out << "*Element,type=C3D8,elset=E\n";
int e = 0;
for (int z = 0; z < 2; ++z)
for (int y = 0; y < 2; ++y)
for (int x = 0; x < 2; ++x)
out << ++e << ',' << id(x, y, z) << ',' << id(x + 1, y, z) << ','
<< id(x + 1, y + 1, z) << ',' << id(x, y + 1, z) << ','
<< id(x, y, z + 1) << ',' << id(x + 1, y, z + 1) << ','
<< id(x + 1, y + 1, z + 1) << ',' << id(x, y + 1, z + 1) << '\n';
out << "*Solid Section,elset=E,material=M\n,\n*End Part\n"
"*Assembly,name=A\n*Instance,name=I,part=P\n*End Instance\n"
"*End Assembly\n*Material,name=M\n*Elastic\n120,.25\n"
"*Step,name=Patch\n*Static\n1,1,1,1\n*Boundary\n";
for (int z = 0; z < 3; ++z)
for (int y = 0; y < 3; ++y)
for (int x = 0; x < 3; ++x) {
if (x == 1 && y == 1 && z == 1) continue;
out << id(x, y, z) << ",1,1," << 0.02 * x + 0.03 * y << '\n';
out << id(x, y, z) << ",2,3,0\n";
}
out << "*End Step\n";
auto domain = Map(out.str());
ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
const auto& frame = result.Value().state.CommittedState();
ASSERT_EQ(frame.solid_results.size(), 64U);
EXPECT_NEAR(frame.displacement[3 * 13], 0.05, 1e-13);
EXPECT_NEAR(frame.displacement[3 * 13 + 1], 0.0, 1e-13);
EXPECT_NEAR(frame.displacement[3 * 13 + 2], 0.0, 1e-13);
EXPECT_NEAR(frame.physical_strain_energy, 8 * 0.0504, 1e-12);
for (std::size_t g = 0; g < 64; ++g) {
EXPECT_EQ(frame.solid_results[g].element, g / 8);
EXPECT_EQ(frame.solid_results[g].integration_point,
static_cast<int>(g % 8 + 1));
EXPECT_NEAR(frame.solid_results[g].cauchy_stress[0], 2.88, 1e-12);
EXPECT_NEAR(frame.solid_results[g].cauchy_stress[3], 1.44, 1e-12);
}
fesa::TbbParallelFor parallel;
fesa::MklPardisoSolver linear;
fesa::MklPardisoTangentSolver tangent;
auto parallel_result =
fesa::SolidStaticProcedureExecutor(parallel, linear, tangent)
.Execute(domain.Value());
ASSERT_TRUE(parallel_result.HasValue());
const auto& parallel_frame = parallel_result.Value().state.CommittedState();
EXPECT_DOUBLE_EQ(parallel_frame.physical_strain_energy,
frame.physical_strain_energy);
for (std::size_t i = 0; i < frame.displacement.Size(); ++i) {
EXPECT_DOUBLE_EQ(parallel_frame.displacement[i], frame.displacement[i]);
EXPECT_DOUBLE_EQ(parallel_frame.internal_force[i], frame.internal_force[i]);
EXPECT_DOUBLE_EQ(parallel_frame.residual[i], frame.residual[i]);
}
for (std::size_t g = 0; g < 64; ++g) {
EXPECT_EQ(parallel_frame.solid_results[g].cauchy_stress,
frame.solid_results[g].cauchy_stress);
EXPECT_EQ(parallel_frame.solid_results[g].strain,
frame.solid_results[g].strain);
}
}
TEST(SolidStaticProcedure,
LaterInvertedPrescribedTrialPreservesThreeAcceptedFrames) {
auto domain = Map(Deck(AffineBoundary(-1.2), "", true));
ASSERT_TRUE(domain.HasValue());
fesa::SerialParallelFor serial;
fesa::MklPardisoSolver linear;
fesa::MklPardisoTangentSolver tangent;
fesa::SolidStaticState state;
auto status = fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
.ExecuteInto(domain.Value(), state);
ASSERT_FALSE(status.IsOk());
ASSERT_EQ(state.AcceptedFrames().size(), 3U);
EXPECT_DOUBLE_EQ(state.CommittedState().load_factor, 0.75);
EXPECT_NEAR(state.CommittedState().displacement[3], -0.9, 1e-15);
ASSERT_EQ(state.CommittedState().solid_results.size(), 8U);
EXPECT_EQ(state.IterationRecords().back().increment_index, 4U);
EXPECT_FALSE(state.IterationRecords().back().accepted);
}
TEST(SolidStaticProcedure, PartiallySupportedRigidMechanismAlsoFailsUnloaded) {
for (const auto& load : {std::string{}, Tension(1.2)}) {
auto domain = Map(Deck("1,1,3,0\n", load));
ASSERT_TRUE(domain.HasValue());
auto result = ExecuteSolid(domain.Value());
EXPECT_FALSE(result.HasValue());
if (!result.HasValue())
EXPECT_EQ(result.GetStatus().Category(), fesa::FailureCategory::kSolver);
}
}
TEST(SolidStaticProcedure, SerialAndTbbFramesAndFailureEvidenceAreIdentical) {
auto domain = Map(Deck(Supports(), Tension(1.2), true));
ASSERT_TRUE(domain.HasValue());
fesa::SerialParallelFor serial;
fesa::TbbParallelFor parallel;
fesa::MklPardisoSolver linear_a, linear_b;
fesa::MklPardisoTangentSolver tangent_a, tangent_b;
auto a = fesa::SolidStaticProcedureExecutor(serial, linear_a, tangent_a)
.Execute(domain.Value());
auto b = fesa::SolidStaticProcedureExecutor(parallel, linear_b, tangent_b)
.Execute(domain.Value());
ASSERT_TRUE(a.HasValue());
ASSERT_TRUE(b.HasValue());
ASSERT_EQ(a.Value().state.AcceptedFrames().size(),
b.Value().state.AcceptedFrames().size());
for (std::size_t n = 0; n < a.Value().state.AcceptedFrames().size(); ++n) {
const auto& af = a.Value().state.AcceptedFrames()[n];
const auto& bf = b.Value().state.AcceptedFrames()[n];
EXPECT_DOUBLE_EQ(af.physical_strain_energy, bf.physical_strain_energy);
EXPECT_EQ(af.global_equilibrium, bf.global_equilibrium);
for (std::size_t i = 0; i < 24; ++i) {
EXPECT_DOUBLE_EQ(af.displacement[i], bf.displacement[i]);
EXPECT_DOUBLE_EQ(af.residual[i], bf.residual[i]);
EXPECT_DOUBLE_EQ(af.reaction[i], bf.reaction[i]);
}
for (std::size_t g = 0; g < 8; ++g) {
EXPECT_EQ(af.solid_results[g].cauchy_stress,
bf.solid_results[g].cauchy_stress);
EXPECT_EQ(af.solid_results[g].strain, bf.solid_results[g].strain);
}
}
auto invalid = Map(Deck(AffineBoundary(-1.2), "", true));
ASSERT_TRUE(invalid.HasValue());
fesa::SolidStaticState sa, sb;
auto as = fesa::SolidStaticProcedureExecutor(serial, linear_a, tangent_a)
.ExecuteInto(invalid.Value(), sa);
auto bs = fesa::SolidStaticProcedureExecutor(parallel, linear_b, tangent_b)
.ExecuteInto(invalid.Value(), sb);
ASSERT_FALSE(as.IsOk());
ASSERT_FALSE(bs.IsOk());
EXPECT_EQ(as.Diagnostics()[0].code, bs.Diagnostics()[0].code);
EXPECT_EQ(sa.AcceptedFrames().size(), sb.AcceptedFrames().size());
EXPECT_EQ(sa.IterationRecords().back().stage,
sb.IterationRecords().back().stage);
EXPECT_DOUBLE_EQ(sa.CommittedState().physical_strain_energy,
sb.CommittedState().physical_strain_energy);
}
TEST(SolidStaticProcedure,
InvalidRecoveryAndRepeatedExecutionKeepAcceptedFrame) {
auto domain = Map(Deck(Supports(), Tension(1.2), true));
ASSERT_TRUE(domain.HasValue());
fesa::SerialParallelFor serial;
fesa::MklPardisoSolver linear;
FailedTangent tangent{FailedTangent::Mode::kBackend};
fesa::SolidStaticState state;
ASSERT_FALSE(fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
.ExecuteInto(domain.Value(), state)
.IsOk());
ASSERT_EQ(state.AcceptedFrames().size(), 1U);
const auto valid = state.CommittedState();
for (int mutation = 0; mutation < 5; ++mutation) {
auto bad = valid;
// Supply the next valid identity so rejection must examine recovery.
bad.frame_index = 1U;
bad.increment_index = 2U;
bad.load_factor = .5;
if (mutation == 0)
bad.solid_results[0].cauchy_stress[2] =
std::numeric_limits<double>::infinity();
if (mutation == 1) bad.solid_results[0].integration_point = 2;
if (mutation == 2) bad.solid_results[0].reference_coordinates[0] += 1;
if (mutation == 3) bad.reaction[0] += 1;
if (mutation == 4) bad.nodal_source_ids[0].source_label = 999;
EXPECT_FALSE(state.CommitAcceptedFrame(bad).IsOk());
ASSERT_EQ(state.AcceptedFrames().size(), 1U);
EXPECT_DOUBLE_EQ(state.CommittedState().physical_strain_energy,
valid.physical_strain_energy);
EXPECT_EQ(state.CommittedState().solid_results[0].cauchy_stress,
valid.solid_results[0].cauchy_stress);
}
EXPECT_FALSE(fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
.ExecuteInto(domain.Value(), state)
.IsOk());
ASSERT_EQ(state.AcceptedFrames().size(), 1U);
}
TEST(SolidStaticProcedure, StressedIndefiniteFreeTangentUsesCapableBackend) {
// F=diag(.2,1,1): the homogeneous transverse right-face mode has
// second variation mu + S11 = 48 - 69.12 < 0. It is not an SPD tangent.
std::ostringstream boundary;
std::ostringstream load;
for (int a = 0; a < 8; ++a) {
boundary << a + 1 << ",1,1," << -.8 * kX[a][0] << '\n';
boundary << a + 1 << ",3,3,0\n";
if (kX[a][0] == 0) boundary << a + 1 << ",2,2,0\n";
// P22 = S22 = lambda*E11 = -23.04; each Y-face node gets P22/4.
// Hold the transverse affine state with its independently known traction.
else
load << a + 1 << ",2," << (kX[a][1] == 0 ? 5.76 : -5.76) << '\n';
}
auto domain = Map(Deck(boundary.str(), load.str(), true, 1.0));
ASSERT_TRUE(domain.HasValue());
fesa::SerialParallelFor serial;
fesa::MklPardisoSolver linear;
ObservedIndefiniteTangent tangent;
auto result = fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
.Execute(domain.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
EXPECT_NEAR(tangent.rayleigh_sum, -21.12, 1e-12);
ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 1U);
const auto& frame = result.Value().state.CommittedState();
for (std::size_t a = 0; a < 8; ++a)
EXPECT_NEAR(frame.displacement[3 * a + 1], 0.0, 1e-13);
EXPECT_NEAR(frame.solid_results[0].strain[0], -.48, 1e-13);
}
TEST(SolidStaticProcedure,
FullyConstrainedZeroDisplacementRetainsSupportLoadResidual) {
auto domain = Map(Deck(AffineBoundary(0.0), "1,1,7\n"));
ASSERT_TRUE(domain.HasValue());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue());
const auto& frame = result.Value().state.CommittedState();
EXPECT_DOUBLE_EQ(frame.internal_force[0], 0.0);
EXPECT_DOUBLE_EQ(frame.residual[0], -7.0);
EXPECT_DOUBLE_EQ(frame.reaction[0], -7.0);
for (double value : frame.global_equilibrium) EXPECT_DOUBLE_EQ(value, 0.0);
EXPECT_DOUBLE_EQ(frame.physical_strain_energy, 0.0);
ASSERT_EQ(frame.solid_results.size(), 8U);
}
TEST(SolidStaticProcedure, UniformPrescribedTranslationHasExactlyZeroRecovery) {
// A rigid translation has no strain, stress, energy, internal force or RF.
// Binary-exact offsets independently expose cancellation in global K*u.
const double translation[3] = {2.75, -3.25, 1.5};
std::ostringstream boundary;
for (int node = 1; node <= 8; ++node)
for (int component = 0; component < 3; ++component)
boundary << node << ',' << component + 1 << ',' << component + 1 << ','
<< translation[component] << '\n';
auto domain = Map(Deck(boundary.str()));
ASSERT_TRUE(domain.HasValue());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
const auto& frame = result.Value().state.CommittedState();
ASSERT_EQ(frame.displacement.Size(), 24U);
for (std::size_t i = 0; i < 24U; ++i) {
EXPECT_DOUBLE_EQ(frame.displacement[i], translation[i % 3U]);
EXPECT_DOUBLE_EQ(frame.internal_force[i], 0.0);
EXPECT_DOUBLE_EQ(frame.residual[i], 0.0);
EXPECT_DOUBLE_EQ(frame.reaction[i], 0.0);
}
EXPECT_DOUBLE_EQ(frame.physical_strain_energy, 0.0);
ASSERT_EQ(frame.solid_results.size(), 8U);
for (const auto& point : frame.solid_results) {
for (double strain : point.strain) EXPECT_DOUBLE_EQ(strain, 0.0);
for (double stress : point.cauchy_stress) EXPECT_DOUBLE_EQ(stress, 0.0);
}
}
TEST(SolidStaticProcedure, MixedStrainAndTranslationPreserveAnalyticRecovery) {
const double translation[3] = {2.75, -3.25, 1.5};
// Binary-exact affine strain: ex=1/32, engineering shear=1/16.
// Exact constant stress is [4.5,1.5,1.5,3,0,0] for E=120, nu=.25.
std::ostringstream translated_boundary;
translated_boundary.precision(17);
for (int a = 0; a < 8; ++a) {
const double u[3] = {translation[0] + kX[a][0] / 32 + kX[a][1] / 16,
translation[1], translation[2]};
for (int component = 0; component < 3; ++component)
translated_boundary << a + 1 << ',' << component + 1 << ','
<< component + 1 << ',' << u[component] << '\n';
}
auto original_domain = Map(Deck(AffineBoundary(1.0 / 32, 1.0 / 16)));
auto translated_domain = Map(Deck(translated_boundary.str()));
ASSERT_TRUE(original_domain.HasValue());
ASSERT_TRUE(translated_domain.HasValue());
auto original = ExecuteSolid(original_domain.Value());
auto translated = ExecuteSolid(translated_domain.Value());
ASSERT_TRUE(original.HasValue()) << Message(original.GetStatus());
ASSERT_TRUE(translated.HasValue()) << Message(translated.GetStatus());
const auto& a = original.Value().state.CommittedState();
const auto& b = translated.Value().state.CommittedState();
const double stress[6] = {4.5, 1.5, 1.5, 3.0, 0.0, 0.0};
for (const auto* frame : {&a, &b}) {
EXPECT_NEAR(frame->physical_strain_energy, 0.1640625, 1e-13);
for (const auto& point : frame->solid_results) {
EXPECT_NEAR(point.strain[0], 1.0 / 32, 1e-15);
EXPECT_NEAR(point.strain[3], 1.0 / 32, 1e-15);
for (std::size_t c = 0; c < 6U; ++c)
EXPECT_NEAR(point.cauchy_stress[c], stress[c], 1e-13);
}
for (std::size_t node = 0; node < 8U; ++node) {
// Independent constant-stress surface traction integrals: each incident
// unit face contributes its outward normal times sigma/4 at this node.
const double nx = (2 * kX[node][0] - 1) / 4;
const double ny = (2 * kX[node][1] - 1) / 4;
const double nz = (2 * kX[node][2] - 1) / 4;
const double exact_force[3] = {4.5 * nx + 3 * ny, 3 * nx + 1.5 * ny,
1.5 * nz};
for (std::size_t c = 0; c < 3U; ++c) {
EXPECT_NEAR(frame->internal_force[3 * node + c], exact_force[c], 1e-13);
EXPECT_NEAR(frame->reaction[3 * node + c], exact_force[c], 1e-13);
}
}
}
EXPECT_DOUBLE_EQ(a.physical_strain_energy, b.physical_strain_energy);
for (std::size_t i = 0; i < 24U; ++i) {
EXPECT_DOUBLE_EQ(a.internal_force[i], b.internal_force[i]);
EXPECT_DOUBLE_EQ(a.reaction[i], b.reaction[i]);
EXPECT_DOUBLE_EQ(a.residual[i], b.residual[i]);
}
for (std::size_t g = 0; g < 8U; ++g) {
EXPECT_EQ(a.solid_results[g].strain, b.solid_results[g].strain);
EXPECT_EQ(a.solid_results[g].cauchy_stress,
b.solid_results[g].cauchy_stress);
}
}
TEST(SolidStaticProcedure, PureLargeRigidRotationEndsAtZeroStrainEnergy) {
auto domain = Map(Deck(AffineBoundary(0.0, 0.0, 1.2), "", true));
ASSERT_TRUE(domain.HasValue());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
const auto& frame = result.Value().state.CommittedState();
EXPECT_NEAR(frame.physical_strain_energy, 0.0, 1e-25);
for (const auto& point : frame.solid_results) {
for (double value : point.strain) EXPECT_NEAR(value, 0.0, 1e-14);
for (double value : point.cauchy_stress) EXPECT_NEAR(value, 0.0, 1e-12);
}
}
TEST(SolidStaticProcedure, OrdinaryStaticKeepsExistingFinalShortFactorPath) {
std::string deck = Deck(AffineBoundary(.03), "", true);
const std::string from = "*Static,direct\n0.25,1";
const auto position = deck.find(from);
ASSERT_NE(position, std::string::npos);
deck.replace(position, from.size(), "*Static\n.3,1,.1,.3");
auto domain = Map(deck);
ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 4U);
const double factors[4] = {.3, .6, .9, 1.0};
for (std::size_t i = 0; i < 4; ++i)
EXPECT_NEAR(result.Value().state.AcceptedFrames()[i].load_factor,
factors[i], 1e-15);
EXPECT_DOUBLE_EQ(result.Value().state.CommittedState().load_factor, 1.0);
}
namespace {
// Actual slender 3D cantilever. The beam limit is a useful trend oracle,
// not an assertion that the coarse 3D discretization equals Euler theory.
std::string Cantilever(int nx, double poisson) {
const int ny = 2, nz = 2;
const auto id = [nx](int x, int y, int z) {
return 1 + x + (nx + 1) * (y + 3 * z);
};
std::ostringstream out;
out << "*Part,name=P\n*Node\n";
for (int z = 0; z <= nz; ++z)
for (int y = 0; y <= ny; ++y)
for (int x = 0; x <= nx; ++x)
out << id(x, y, z) << ',' << 8.0 * x / nx << ',' << .5 * y << ','
<< .5 * z << '\n';
out << "*Element,type=C3D8,elset=E\n";
int e = 0;
for (int z = 0; z < nz; ++z)
for (int y = 0; y < ny; ++y)
for (int x = 0; x < nx; ++x)
out << ++e << ',' << id(x, y, z) << ',' << id(x + 1, y, z) << ','
<< id(x + 1, y + 1, z) << ',' << id(x, y + 1, z) << ','
<< id(x, y, z + 1) << ',' << id(x + 1, y, z + 1) << ','
<< id(x + 1, y + 1, z + 1) << ',' << id(x, y + 1, z + 1) << '\n';
out << "*Solid Section,elset=E,material=M\n,\n*End Part\n"
"*Assembly,name=A\n*Instance,name=I,part=P\n*End Instance\n"
"*End Assembly\n*Material,name=M\n*Elastic\n120,"
<< poisson << '\n';
out << "*Step,name=Bending\n*Static\n1,1,1,1\n*Boundary\n";
for (int z = 0; z <= nz; ++z)
for (int y = 0; y <= ny; ++y) out << id(0, y, z) << ",1,3,0\n";
out << "*Cload\n";
for (int z = 0; z <= nz; ++z)
for (int y = 0; y <= ny; ++y) {
const double weight = (y == 1 ? .5 : .25) * (z == 1 ? .5 : .25);
out << id(nx, y, z) << ",2," << 1e-4 * weight << '\n';
}
out << "*End Step\n";
return out.str();
}
} // namespace
TEST(SolidStaticProcedure,
BoundedSpatialRefinementReportsFullIntegrationLocking) {
const double beam_limit = 4 * 1e-4 * 8 * 8 * 8 / 120;
double previous[2] = {0, 0};
double final[2] = {0, 0};
for (int material = 0; material < 2; ++material) {
const double poisson = material == 0 ? .25 : .49;
for (int nx : {2, 4, 8}) {
auto domain = Map(Cantilever(nx, poisson));
ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
const auto& frame = result.Value().state.CommittedState();
// End-center index x=nx,y=1,z=1 in node-major source ordering.
const auto node = static_cast<std::size_t>(nx + 4 * (nx + 1));
const double tip = frame.displacement[3 * node + 1];
std::cout << "solid-static-refinement nu=" << poisson << " nx=" << nx
<< " tip=" << tip << " beam-limit=" << beam_limit << '\n';
EXPECT_GT(tip, previous[material]);
EXPECT_LT(tip, 1.05 * beam_limit);
previous[material] = tip;
final[material] = tip;
}
}
// Near-incompressibility remains stiffer with this full-integration mesh.
EXPECT_LT(final[1], final[0]);
EXPECT_GT(std::abs(final[1] - beam_limit), std::abs(final[0] - beam_limit));
}
TEST(SolidStaticProcedure,
LaterIterationNonfiniteAndGeometryFailureKeepEveryAcceptedQuantity) {
auto domain = Map(Deck(Supports(), Tension(1.2), true));
ASSERT_TRUE(domain.HasValue());
auto successful = ExecuteSolid(domain.Value());
ASSERT_TRUE(successful.HasValue());
const auto& expected = successful.Value().state.AcceptedFrames()[0];
ASSERT_EQ(expected.iteration_count, 2U);
for (auto mode :
{FailedTangent::Mode::kGeometry, FailedTangent::Mode::kIteration,
FailedTangent::Mode::kNonfinite}) {
fesa::SerialParallelFor serial;
fesa::MklPardisoSolver linear;
FailedTangent tangent{mode, true};
fesa::SolidStaticState state;
auto status = fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
.ExecuteInto(domain.Value(), state);
ASSERT_FALSE(status.IsOk());
ASSERT_EQ(state.AcceptedFrames().size(), 1U);
const auto& actual = state.CommittedState();
EXPECT_EQ(actual.global_equilibrium, expected.global_equilibrium);
EXPECT_DOUBLE_EQ(actual.physical_strain_energy,
expected.physical_strain_energy);
for (std::size_t i = 0; i < 24; ++i) {
EXPECT_DOUBLE_EQ(actual.displacement[i], expected.displacement[i]);
EXPECT_DOUBLE_EQ(actual.internal_force[i], expected.internal_force[i]);
EXPECT_DOUBLE_EQ(actual.external_force[i], expected.external_force[i]);
EXPECT_DOUBLE_EQ(actual.residual[i], expected.residual[i]);
EXPECT_DOUBLE_EQ(actual.reaction[i], expected.reaction[i]);
}
for (std::size_t g = 0; g < 8; ++g) {
EXPECT_EQ(actual.solid_results[g].strain,
expected.solid_results[g].strain);
EXPECT_EQ(actual.solid_results[g].cauchy_stress,
expected.solid_results[g].cauchy_stress);
}
EXPECT_EQ(state.IterationRecords().back().increment_index, 2U);
EXPECT_FALSE(state.IterationRecords().back().accepted);
}
}
TEST(SolidStaticProcedure,
ExecutesBothUnmodifiedStaticSourceDecksThroughDirectApi) {
auto root = std::filesystem::path{__FILE__}.parent_path();
for (int i = 0; i < 3; ++i) root = root.parent_path();
for (const std::string name :
{"hexa_linearstatic.inp", "hexa_nonlinearstatic.inp"}) {
auto syntax =
fesa::AbaqusInputReader{}.Read(root / "reference/hexa" / name);
ASSERT_TRUE(syntax.HasValue()) << Message(syntax.GetStatus());
auto domain = fesa::AbaqusDomainMapper{}.Map(syntax.Value());
ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
auto result = ExecuteSolid(domain.Value());
ASSERT_TRUE(result.HasValue())
<< name << ' ' << Message(result.GetStatus());
const auto& state = result.Value().state;
const auto& frame = state.CommittedState();
EXPECT_DOUBLE_EQ(frame.load_factor, 1.0);
EXPECT_EQ(frame.increment_index, 10U);
EXPECT_EQ(frame.displacement.Size(), 442U * 3U);
EXPECT_EQ(frame.solid_results.size(), 300U * 8U);
EXPECT_EQ(state.AcceptedFrames().size(),
domain.Value().Steps()[0].IsNonlinearStatic() ? 10U : 1U);
std::cout << "solid-static-source " << name
<< " frames=" << state.AcceptedFrames().size()
<< " iterations=" << state.IterationRecords().size()
<< " energy=" << frame.physical_strain_energy << '\n';
}
}
tests/unit/analysis/solid_static_state_test.cpp
+60
View File
@@ -0,0 +1,60 @@
#include "fesa/analysis/solid_static_state.h"
#include <gtest/gtest.h>
#include <limits>
TEST(SolidStaticControls, PurePrescribedUsesAbsoluteForceAndOnlyTwoChannels) {
fesa::Vector residual{1U};
fesa::Vector internal{1U};
fesa::Vector external{1U};
fesa::Vector correction{1U};
fesa::Vector increment{1U};
residual[0] = internal[0] = 0.0009;
correction[0] = 0.0009;
auto result = fesa::EvaluateSolidStaticConvergence(
residual, internal, external, correction, increment);
ASSERT_TRUE(result.HasValue());
EXPECT_EQ(result.Value().force_residual.branch,
fesa::NewtonConvergenceBranch::kZeroScale);
EXPECT_DOUBLE_EQ(result.Value().force_residual.denominator_scale, 0.0);
EXPECT_EQ(result.Value().force_residual.threshold_unit,
fesa::NewtonChannelUnit::kForce);
EXPECT_EQ(result.Value().translation_correction.threshold_unit,
fesa::NewtonChannelUnit::kLength);
EXPECT_TRUE(result.Value().force_residual.passed);
EXPECT_TRUE(result.Value().translation_correction.passed);
internal[0] = residual[0] = 0.0011;
EXPECT_FALSE(fesa::EvaluateSolidStaticConvergence(
residual, internal, external, correction, increment)
.Value()
.force_residual.passed);
}
TEST(SolidStaticControls, RelativeLimitsRawNormsAndNonfiniteFailure) {
fesa::Vector r{1U}, f{1U}, p{1U}, du{1U}, u{1U};
r[0] = 1.0;
f[0] = 100.0;
p[0] = 99.0;
du[0] = 0.1;
u[0] = 10.0;
auto result = fesa::EvaluateSolidStaticConvergence(r, f, p, du, u);
ASSERT_TRUE(result.HasValue());
EXPECT_DOUBLE_EQ(result.Value().force_residual.numerator_norm, 1.0);
EXPECT_DOUBLE_EQ(result.Value().force_residual.denominator_scale, 100.0);
EXPECT_EQ(result.Value().force_residual.branch,
fesa::NewtonConvergenceBranch::kRelative);
EXPECT_EQ(result.Value().force_residual.threshold_unit,
fesa::NewtonChannelUnit::kDimensionless);
EXPECT_TRUE(result.Value().force_residual.passed);
EXPECT_TRUE(result.Value().translation_correction.passed);
r[0] = 1.0001;
EXPECT_FALSE(fesa::EvaluateSolidStaticConvergence(r, f, p, du, u)
.Value()
.force_residual.passed);
p[0] = std::numeric_limits<double>::infinity();
EXPECT_FALSE(fesa::EvaluateSolidStaticConvergence(r, f, p, du, u).HasValue());
EXPECT_FALSE(
fesa::EvaluateSolidStaticConvergence(r, f, fesa::Vector{2U}, du, u)
.HasValue());
}