diff --git a/include/fesa/analysis/solid_static_procedure.h b/include/fesa/analysis/solid_static_procedure.h new file mode 100644 index 0000000..d36b418 --- /dev/null +++ b/include/fesa/analysis/solid_static_procedure.h @@ -0,0 +1,39 @@ +#ifndef FESA_ANALYSIS_SOLID_STATIC_PROCEDURE_H_ +#define FESA_ANALYSIS_SOLID_STATIC_PROCEDURE_H_ + +#include "fesa/analysis/solid_static_state.h" + +namespace fesa { +class LinearSolver; +class ParallelFor; +class TangentSolver; + +/// @brief Owns successful static frames and sorted model diagnostics. +struct SolidStaticExecutionResult { + SolidStaticState state; + std::vector diagnostics; +}; + +/// @brief Executes linear K0 or fixed-factor TL StVK static analysis. +class SolidStaticProcedureExecutor { + public: + /// @brief Uses non-owned adapters, which must outlive this executor. + /// @param tangent_solver Must support indefinite stressed tangents. + SolidStaticProcedureExecutor(const ParallelFor& parallel_for, + LinearSolver& linear_solver, + TangentSolver& tangent_solver); + /// @brief Returns validated results; no I/O or final-file replacement. + Result Execute(const Domain& domain) const; + /// @brief Executes into a fresh observer retaining accepted history on + /// failure. + /// @note A previously initialized observer is rejected without mutation. + Status ExecuteInto(const Domain& domain, SolidStaticState& state) const; + + private: + const ParallelFor& parallel_for_; + LinearSolver& linear_solver_; + TangentSolver& tangent_solver_; +}; +} // namespace fesa + +#endif // FESA_ANALYSIS_SOLID_STATIC_PROCEDURE_H_ diff --git a/include/fesa/analysis/solid_static_state.h b/include/fesa/analysis/solid_static_state.h new file mode 100644 index 0000000..d45eae2 --- /dev/null +++ b/include/fesa/analysis/solid_static_state.h @@ -0,0 +1,115 @@ +#ifndef FESA_ANALYSIS_SOLID_STATIC_STATE_H_ +#define FESA_ANALYSIS_SOLID_STATIC_STATE_H_ + +#include +#include +#include +#include + +#include "fesa/analysis/nonlinear_static_control.h" +#include "fesa/results/solid_result_records.h" + +namespace fesa { + +/// @brief Keeps raw current-arm balance and its free-residual explanation. +/// @note Internal evidence only; no additional persistent result fields. +struct SolidStaticBalanceEvidence { + std::array internal_force_sum{}; + std::array reference_internal_moment{}; + std::array current_internal_moment{}; + std::array free_residual_moment{}; + std::array linearization_moment{}; + std::array internal_force_scale{}; + std::array reference_moment_scale{}; + std::array current_moment_scale{}; +}; + +/// @brief Owns a complete accepted static frame in stable node/source order. +/// @note Vectors are node-major UX/UY/UZ; reaction is constrained residual. +struct SolidStaticFrame { + std::uint32_t frame_index{0U}; + std::uint32_t increment_index{0U}; + double load_factor{0.0}; + std::uint32_t iteration_count{0U}; + bool is_factor_one{false}; + std::vector nodal_internal_order; + std::vector nodal_source_ids; + Vector displacement{0U}; + Vector internal_force{0U}; + Vector external_force{0U}; + Vector residual{0U}; + Vector reaction{0U}; + std::vector solid_results; + double physical_strain_energy{0.0}; + // Applied plus support forces, then moments about global zero at X+u. + std::array global_equilibrium{}; + SolidStaticBalanceEvidence balance; +}; + +/// @brief Records two physical static channels, including failed iterations. +struct SolidStaticIterationRecord { + std::uint32_t increment_index{0U}; + double load_factor{0.0}; + std::uint32_t iteration{0U}; + NewtonChannelDecision force_residual{"residual_force", + 0.0, + 0.0, + NewtonConvergenceBranch::kZeroScale, + kNewtonZeroScaleThreshold, + NewtonChannelUnit::kForce, + false}; + NewtonChannelDecision translation_correction{ + "correction_translation", + 0.0, + 0.0, + NewtonConvergenceBranch::kZeroScale, + kNewtonZeroScaleThreshold, + NewtonChannelUnit::kLength, + false}; + bool convergence_evaluated{false}; + bool accepted{false}; + std::string diagnostic_code; + std::string stage; + std::string failure_reason; +}; + +/// @brief Evaluates translation-only static criteria with input force/length. +/// @note Exact-zero external force selects the absolute residual branch. +Result EvaluateSolidStaticConvergence( + const Vector& residual, const Vector& internal, const Vector& external, + const Vector& correction, const Vector& increment); + +/// @brief Owns committed static state independently of caller trial work. +class SolidStaticState { + public: + /// @brief Creates an uninitialized observer usable with ExecuteInto(). + SolidStaticState() = default; + /// @brief Copies immutable source inventory and the zero-factor baseline. + static Result Create(const Domain& domain, + const DofManager& dofs); + /// @brief Validates complete recovery, then appends and replaces atomically. + /// @note The executor certifies trial geometry before this state boundary. + Status CommitAcceptedFrame(const SolidStaticFrame& frame); + /// @brief Returns authoritative accepted frames; baseline is not a frame. + const std::vector& AcceptedFrames() const noexcept; + /// @brief Returns baseline or the most recently accepted complete state. + const SolidStaticFrame& CommittedState() const noexcept; + /// @brief Returns all iteration evidence, including the final failed trial. + const std::vector& IterationRecords() + const noexcept; + + private: + friend class SolidStaticProcedureExecutor; + bool initialized_{false}; + bool nonlinear_{false}; + std::vector reference_points_; + std::vector constrained_dofs_; + Vector prescribed_values_{0U}; + SolidStaticFrame committed_; + std::vector accepted_frames_; + std::vector iterations_; +}; + +} // namespace fesa + +#endif // FESA_ANALYSIS_SOLID_STATIC_STATE_H_ diff --git a/src/fesa/analysis/solid_static_procedure.cpp b/src/fesa/analysis/solid_static_procedure.cpp new file mode 100644 index 0000000..a18b7eb --- /dev/null +++ b/src/fesa/analysis/solid_static_procedure.cpp @@ -0,0 +1,471 @@ +#include "fesa/analysis/solid_static_procedure.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/assembly/load_assembler.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/assembly/sparse_assembler.h" +#include "fesa/constraints/essential_constraint_policy.h" +#include "fesa/elements/element_factory.h" +#include "fesa/elements/hex8_element.h" +#include "fesa/math/sparse_matrix.h" +#include "fesa/solvers/linear/linear_solver.h" +#include "fesa/solvers/linear/tangent_solver.h" + +namespace fesa { +namespace { +constexpr std::uint32_t kMaximumIterations = 100U; +constexpr double kRoundoffBudget = + 512.0 * std::numeric_limits::epsilon(); + +struct RuntimeSolid { + const Hex8Element* element; + std::vector scatter; +}; +struct SolidAssembly { + Vector force; + SparseMatrix tangent; + std::vector points; + double energy; +}; + +/// @brief Attaches immutable source-step context to a procedure failure. +Status Failure(const Domain& domain, const char* code, const char* message, + FailureCategory category = FailureCategory::kSolver) { + return Status::Failure(category, {{Severity::kError, + code, + {domain.SourcePath(), 0U}, + "STATIC", + domain.SourceContentIdentity(), + message}}); +} +/// @brief Preserves kernel source/entity evidence for failed trial recovery. +Status TrialFailure(const Status& status) { + return Status::Failure(FailureCategory::kSolver, status.Diagnostics()); +} +/// @brief Returns component products with uncancelled arithmetic scales. +std::array Cross(const std::array& a, + const std::array& b, + std::array* scale = nullptr) { + std::array result{}; + for (std::size_t i = 0U; i < 3U; ++i) { + const auto j = (i + 1U) % 3U; + const auto k = (i + 2U) % 3U; + result[i] = a[j] * b[k] - a[k] * b[j]; + if (scale != nullptr) + (*scale)[i] += std::abs(a[j] * b[k]) + std::abs(a[k] * b[j]); + } + return result; +} +/// @brief Tests a zero identity at its uncancelled same-unit operation scale. +bool RoundoffZero(double value, double scale) { + return std::isfinite(value) && std::isfinite(scale) && + (scale == 0.0 ? value == 0.0 + : std::abs(value) <= kRoundoffBudget * scale); +} +/// @brief Derives source DIRECT factors or the existing ordinary STATIC path. +Result> Factors(const Domain& domain) { + const auto& step = domain.Steps()[0]; + if (!step.IsDirectStatic()) return BuildFixedLoadFactorSequence(step); + const auto count = step.IncrementCount(); + if (count == 0U || count > step.MaximumIncrements() || + !std::isfinite(step.InitialIncrement()) || + !(step.InitialIncrement() > 0.0) || !std::isfinite(step.TimePeriod()) || + !(step.TimePeriod() > 0.0) || + step.InitialIncrement() * static_cast(count) != step.TimePeriod()) + return Result>::Failure(Failure( + domain, "invalid-solid-static-schedule", + "DIRECT static requires a finite positive integer fixed schedule.", + FailureCategory::kInput)); + std::vector factors; + factors.reserve(count); + double previous = 0.0; + for (std::uint32_t n = 1U; n <= count; ++n) { + const double factor = n == count + ? 1.0 + : static_cast(n) * + step.InitialIncrement() / step.TimePeriod(); + if (!std::isfinite(factor) || !(factor > previous) || factor > 1.0) + return Result>::Failure( + Failure(domain, "invalid-solid-static-schedule", + "DIRECT static factors must strictly increase to one.", + FailureCategory::kInput)); + factors.push_back({n, factor, n == count}); + previous = factor; + if (n == count) break; // Avoid unsigned wrap at the source limit. + } + return Result>::Success(std::move(factors)); +} +/// @brief Builds index-owned trials, then reduces all contributions stably. +Result Assemble(const std::vector& solids, + const DofManager& dofs, + const Vector& displacement, bool nonlinear, + const ParallelFor& parallel) { + std::vector>> responses(solids.size()); + parallel.Execute(solids.size(), [&](std::size_t e) { + Vector local{24U}; + for (std::size_t i = 0U; i < 24U; ++i) + local[i] = displacement[solids[e].scatter[i]]; + responses[e] = nonlinear ? solids[e].element->EvaluateNonlinear(local) + : solids[e].element->EvaluateLinear(local); + }); + Vector force{dofs.FullDofCount()}; + std::vector contributions; + contributions.reserve(solids.size() * 24U * 24U); + std::vector points; + points.reserve(solids.size() * 8U); + double energy = 0.0; + for (std::size_t e = 0U; e < solids.size(); ++e) { + if (!responses[e]->HasValue()) + return Result::Failure( + TrialFailure(responses[e]->GetStatus())); + const auto& response = responses[e]->Value(); + for (std::size_t i = 0U; i < 24U; ++i) { + force[solids[e].scatter[i]] += response.internal_force[i]; + for (std::size_t j = 0U; j < 24U; ++j) + contributions.push_back({solids[e].scatter[i], solids[e].scatter[j], + response.tangent(i, j), e, 24U * i + j}); + } + energy += response.strain_energy; + points.insert(points.end(), response.points.begin(), response.points.end()); + } + auto tangent = + SparseMatrix::FromCoo(dofs.FullDofCount(), dofs.FullDofCount(), + std::move(contributions), dofs.GetSparsePattern()); + if (!tangent.HasValue()) + return Result::Failure(tangent.GetStatus()); + return Result::Success({std::move(force), + std::move(tangent.Value()), + std::move(points), energy}); +} +/// @brief Computes physical residual, projection and both moment-arm +/// identities. +Status CompleteFrame(const Domain& domain, const DofManager& dofs, + bool nonlinear, const SparseMatrix& stiffness, + SolidStaticFrame& frame) { + const std::size_t size = dofs.FullDofCount(); + std::vector constrained(size, false); + for (auto i : dofs.ConstrainedDofs()) constrained[i] = true; + frame.residual = Vector{size}; + frame.reaction = Vector{size}; + for (std::size_t i = 0U; i < size; ++i) { + frame.residual[i] = frame.internal_force[i] - frame.external_force[i]; + frame.reaction[i] = constrained[i] ? frame.residual[i] : 0.0; + if (!std::isfinite(frame.displacement[i]) || + !std::isfinite(frame.internal_force[i]) || + !std::isfinite(frame.external_force[i]) || + !std::isfinite(frame.residual[i])) + return Failure(domain, "nonfinite-solid-static-recovery", + "Static full-space quantities must remain finite."); + } + if (!std::isfinite(frame.physical_strain_energy)) + return Failure(domain, "nonfinite-solid-static-recovery", + "Recovered static strain energy must be finite."); + std::array force_scale{}; + std::array reference_moment_scale{}; + std::array current_moment_scale{}; + std::array identity_moment_scale{}; + std::array free_force_sum{}; + std::array balance_scale{}; + Vector arithmetic_scale{size}; + // Linear K0 multiplication can cancel large translations within each row. + // Use those original products, not the small resulting internal force. + if (!nonlinear) { + for (std::size_t i = 0U; i < size; ++i) + for (std::size_t k = stiffness.RowOffsets()[i]; + k < stiffness.RowOffsets()[i + 1U]; ++k) + arithmetic_scale[i] += + std::abs(stiffness.Values()[k] * + frame.displacement[stiffness.ColumnIndices()[k]]); + for (auto i : dofs.FreeDofs()) + if (!RoundoffZero( + frame.residual[i], + arithmetic_scale[i] + std::abs(frame.external_force[i]))) + return Failure( + domain, "solid-linear-residual-failed", + "Linear free residual failed its uncancelled roundoff check."); + } + for (std::size_t a = 0U; a < domain.Nodes().size(); ++a) { + const auto& reference = domain.Nodes()[a].coordinates; + std::array current{}; + std::array u{}; + std::array f{}; + std::array e{}; + std::array applied{}; + for (std::size_t i = 0U; i < 3U; ++i) { + const auto full = 3U * a + i; + u[i] = frame.displacement[full]; + current[i] = reference[i] + u[i]; + f[i] = frame.internal_force[full]; + e[i] = constrained[full] ? 0.0 : frame.residual[full]; + applied[i] = frame.external_force[full] + frame.reaction[full]; + frame.balance.internal_force_sum[i] += f[i]; + free_force_sum[i] += e[i]; + force_scale[i] += nonlinear ? std::abs(f[i]) : arithmetic_scale[full]; + frame.global_equilibrium[i] += applied[i]; + balance_scale[i] += std::abs(frame.external_force[full]) + + std::abs(frame.reaction[full]) + std::abs(f[i]) + + std::abs(e[i]); + } + const auto reference_moment = Cross(reference, f, &reference_moment_scale); + const auto current_moment = Cross(current, f, ¤t_moment_scale); + const auto free_moment = Cross(current, e, &identity_moment_scale); + const auto linear_moment = Cross(u, f, &identity_moment_scale); + std::array applied_scale{}; + const auto applied_moment = Cross(current, applied, &applied_scale); + for (std::size_t i = 0U; i < 3U; ++i) { + frame.balance.reference_internal_moment[i] += reference_moment[i]; + frame.balance.current_internal_moment[i] += current_moment[i]; + frame.balance.free_residual_moment[i] += free_moment[i]; + frame.balance.linearization_moment[i] += linear_moment[i]; + frame.global_equilibrium[3U + i] += applied_moment[i]; + balance_scale[3U + i] += applied_scale[i]; + if (!nonlinear) { + const auto j = (i + 1U) % 3U; + const auto k = (i + 2U) % 3U; + reference_moment_scale[i] += + std::abs(reference[j]) * arithmetic_scale[3U * a + k] + + std::abs(reference[k]) * arithmetic_scale[3U * a + j]; + } + } + } + frame.balance.internal_force_scale = force_scale; + frame.balance.reference_moment_scale = reference_moment_scale; + frame.balance.current_moment_scale = current_moment_scale; + for (std::size_t i = 0U; i < 3U; ++i) { + const double physical_moment = + nonlinear ? frame.balance.current_internal_moment[i] + : frame.balance.reference_internal_moment[i]; + const double moment_scale = + nonlinear ? current_moment_scale[i] : reference_moment_scale[i]; + if (!RoundoffZero(frame.balance.internal_force_sum[i], force_scale[i]) || + !RoundoffZero(physical_moment, moment_scale)) + return Failure( + domain, "solid-static-internal-balance-failed", + "Internal force or appropriate moment-arm identity failed."); + const double force_identity = frame.global_equilibrium[i] - + frame.balance.internal_force_sum[i] + + free_force_sum[i]; + const double moment_identity = frame.global_equilibrium[3U + i] - + frame.balance.current_internal_moment[i] + + frame.balance.free_residual_moment[i]; + if (!RoundoffZero(force_identity, balance_scale[i]) || + !RoundoffZero(moment_identity, balance_scale[3U + i] + + current_moment_scale[i] + + identity_moment_scale[i])) + return Failure( + domain, "solid-static-residual-balance-failed", + "Applied/support balance does not explain the free residual."); + if (!std::isfinite(frame.global_equilibrium[i]) || + !std::isfinite(frame.global_equilibrium[3U + i])) + return Failure(domain, "nonfinite-solid-static-recovery", + "Static current-position balance must remain finite."); + } + return Status::Ok(); +} +} // namespace + +SolidStaticProcedureExecutor::SolidStaticProcedureExecutor( + const ParallelFor& parallel_for, LinearSolver& linear_solver, + TangentSolver& tangent_solver) + : parallel_for_(parallel_for), + linear_solver_(linear_solver), + tangent_solver_(tangent_solver) {} + +Status SolidStaticProcedureExecutor::ExecuteInto( + const Domain& domain, SolidStaticState& state) const { + if (state.initialized_) + return Failure(domain, "solid-static-observer-not-fresh", + "Execution requires a fresh static observer."); + if (!domain.IsSolid() || domain.Frequency() != nullptr || + domain.Dynamic() != nullptr || domain.Steps().Size() != 1U) + return Failure(domain, "invalid-solid-static-procedure", + "Solid static execution requires one Solid STATIC step.", + FailureCategory::kModel); + auto model = AnalysisModel::Create(domain); + if (!model.HasValue()) return model.GetStatus(); + std::vector> elements; + ElementView view; + for (auto e : model.Value().ActiveElements()) { + auto element = ElementFactory{}.Create(domain.Elements()[e], domain); + if (!element.HasValue()) return element.GetStatus(); + view.emplace_back(*element.Value()); + elements.push_back(std::move(element.Value())); + } + DofManager dofs; + const Status mapping = dofs.Build(model.Value(), view); + if (!mapping.IsOk()) return mapping; + auto baseline = SolidStaticState::Create(domain, dofs); + if (!baseline.HasValue()) return baseline.GetStatus(); + state = std::move(baseline.Value()); + std::vector solids; + for (const auto& element : elements) { + const auto* hex = dynamic_cast(element.get()); + if (hex == nullptr) + return Failure(domain, "invalid-solid-static-element", + "Solid STATIC requires full-integration Hex8 kernels.", + FailureCategory::kModel); + auto scatter = dofs.ElementScatter(hex->DofLayout()); + if (!scatter.HasValue()) return scatter.GetStatus(); + solids.push_back({hex, std::move(scatter.Value())}); + } + auto stiffness = SparseAssembler::Assemble(view, dofs, parallel_for_); + if (!stiffness.HasValue()) return stiffness.GetStatus(); + const EssentialConstraintPolicy constraints; + const auto& step = model.Value().Step(); + const bool nonlinear = step.IsNonlinearStatic(); + // No support can remove any rigid mode; fail before a perturbed + // factorization. + if (dofs.ConstrainedDofCount() == 0U) + return Failure(domain, "solid-static-rigid-mechanism", + "An unsupported Solid static free system has rigid modes."); + if (!nonlinear) { + auto partition = constraints.Partition(stiffness.Value(), dofs); + if (!partition.HasValue()) return partition.GetStatus(); + Status status = linear_solver_.Factorize(partition.Value().k_ff); + if (!status.IsOk()) return status; + // Preserve the required factorization-before-load-assembly lifecycle. + auto load = LoadAssembler::AssembleFullNodalLoad(model.Value(), dofs); + if (!load.HasValue()) return load.GetStatus(); + auto rhs = LoadAssembler::EffectiveFreeRhs( + load.Value(), partition.Value().k_fc, dofs.PrescribedValues(), dofs); + if (!rhs.HasValue()) return rhs.GetStatus(); + Vector free{dofs.FreeDofCount()}; + status = linear_solver_.Solve(rhs.Value(), free); + if (!status.IsOk()) return status; + SolidStaticFrame frame = state.committed_; + frame.frame_index = 0U; + frame.increment_index = step.IsDirectStatic() ? step.IncrementCount() : 1U; + frame.load_factor = 1.0; + frame.is_factor_one = true; + frame.displacement = + constraints.ReconstructFull(free, dofs.PrescribedValues(), dofs); + auto recovery = + Assemble(solids, dofs, frame.displacement, false, parallel_for_); + if (!recovery.HasValue()) return recovery.GetStatus(); + frame.internal_force = stiffness.Value().Multiply(frame.displacement); + frame.external_force = std::move(load.Value()); + frame.solid_results = std::move(recovery.Value().points); + frame.physical_strain_energy = recovery.Value().energy; + status = CompleteFrame(domain, dofs, false, stiffness.Value(), frame); + return status.IsOk() ? state.CommitAcceptedFrame(frame) : status; + } + auto factors = Factors(domain); + if (!factors.HasValue()) return factors.GetStatus(); + auto base_load = LoadAssembler::AssembleFullNodalLoad(model.Value(), dofs); + if (!base_load.HasValue()) return base_load.GetStatus(); + for (const auto& increment : factors.Value()) { + const Vector start = state.committed_.displacement; + Vector trial = start; + for (std::size_t c = 0U; c < dofs.ConstrainedDofCount(); ++c) + trial[dofs.ConstrainedDofs()[c]] = + increment.load_factor * dofs.PrescribedValues()[c]; + for (std::uint32_t iteration = 1U; iteration <= kMaximumIterations; + ++iteration) { + SolidStaticIterationRecord record; + record.increment_index = increment.increment_index; + record.load_factor = increment.load_factor; + record.iteration = iteration; + record.stage = "entered"; + const auto fail = [&](const Status& failure) { + record.diagnostic_code = failure.Diagnostics().empty() + ? "solid-static-failure" + : failure.Diagnostics()[0].code; + record.failure_reason = failure.Diagnostics().empty() + ? "Static trial failed." + : failure.Diagnostics()[0].message; + state.iterations_.push_back(record); + return failure; + }; + record.stage = "trial-recovery"; + auto assembly = Assemble(solids, dofs, trial, true, parallel_for_); + if (!assembly.HasValue()) return fail(assembly.GetStatus()); + SolidStaticFrame frame = state.committed_; + frame.balance = {}; + frame.global_equilibrium = {}; + frame.frame_index = + static_cast(state.accepted_frames_.size()); + frame.increment_index = increment.increment_index; + frame.load_factor = increment.load_factor; + frame.iteration_count = iteration; + frame.is_factor_one = increment.is_factor_one; + frame.displacement = trial; + frame.internal_force = std::move(assembly.Value().force); + frame.external_force = base_load.Value(); + frame.external_force.Scale(increment.load_factor); + frame.solid_results = std::move(assembly.Value().points); + frame.physical_strain_energy = assembly.Value().energy; + auto status = CompleteFrame(domain, dofs, true, stiffness.Value(), frame); + if (!status.IsOk()) return fail(status); + record.stage = "tangent-partition"; + auto partition = constraints.Partition(assembly.Value().tangent, dofs); + if (!partition.HasValue()) return fail(partition.GetStatus()); + Vector correction{dofs.FreeDofCount()}; + // Empty free is vacuous; prescribed geometry/recovery above still runs. + if (dofs.FreeDofCount() != 0U) { + record.stage = "tangent-factorization"; + status = tangent_solver_.FactorizeFresh(partition.Value().k_ff); + if (!status.IsOk()) return fail(status); + auto rhs = constraints.GatherFree(frame.residual, dofs); + rhs.Scale(-1.0); + record.stage = "correction-solve"; + status = tangent_solver_.SolveCorrection(rhs, correction); + if (!status.IsOk()) return fail(status); + } + Vector current_increment{dofs.FreeDofCount()}; + for (std::size_t f = 0U; f < dofs.FreeDofCount(); ++f) { + const auto full = dofs.FreeDofs()[f]; + current_increment[f] = trial[full] - start[full]; + } + record.stage = "convergence"; + auto convergence = EvaluateSolidStaticConvergence( + constraints.GatherFree(frame.residual, dofs), + constraints.GatherFree(frame.internal_force, dofs), + constraints.GatherFree(frame.external_force, dofs), correction, + current_increment); + if (!convergence.HasValue()) return fail(convergence.GetStatus()); + record.force_residual = convergence.Value().force_residual; + record.translation_correction = + convergence.Value().translation_correction; + record.convergence_evaluated = true; + if (record.force_residual.passed && + record.translation_correction.passed) { + record.stage = "commit"; + status = state.CommitAcceptedFrame(frame); + if (!status.IsOk()) return fail(status); + record.accepted = true; + state.iterations_.push_back(record); + break; + } + if (iteration == kMaximumIterations) + return fail(Failure(domain, "solid-static-nonconvergence", + "Solid STATIC exceeded 100 Newton iterations.")); + for (std::size_t f = 0U; f < dofs.FreeDofCount(); ++f) + trial[dofs.FreeDofs()[f]] += correction[f]; + record.stage = "continue"; + state.iterations_.push_back(record); + } + } + return Status::Ok(); +} +Result SolidStaticProcedureExecutor::Execute( + const Domain& domain) const { + SolidStaticState state; + const Status status = ExecuteInto(domain, state); + if (!status.IsOk()) + return Result::Failure(status); + auto diagnostics = domain.Warnings(); + SortDiagnostics(diagnostics); + return Result::Success( + {std::move(state), std::move(diagnostics)}); +} +} // namespace fesa diff --git a/src/fesa/analysis/solid_static_state.cpp b/src/fesa/analysis/solid_static_state.cpp new file mode 100644 index 0000000..02bfb52 --- /dev/null +++ b/src/fesa/analysis/solid_static_state.cpp @@ -0,0 +1,216 @@ +#include "fesa/analysis/solid_static_state.h" + +#include +#include +#include + +#include "fesa/elements/element_factory.h" +#include "fesa/elements/hex8_element.h" + +namespace fesa { +namespace { +/// @brief Keeps state/control failures separate from input syntax failures. +Status StateFailure(const char* code, const char* message) { + return Status::Failure( + FailureCategory::kSolver, + {{Severity::kError, code, {}, "SOLID_STATIC", {}, message}}); +} +/// @brief Validates all components without a zero clamp or scale floor. +template +bool Finite(const Range& values) { + for (double value : values) + if (!std::isfinite(value)) return false; + return true; +} +/// @brief Validates owned vectors including their physical full-space size. +bool FiniteVector(const Vector& values, std::size_t size) { + if (values.Size() != size) return false; + for (std::size_t i = 0U; i < size; ++i) + if (!std::isfinite(values[i])) return false; + return true; +} +/// @brief Calculates the infinity norm after separate finite validation. +double InfinityNorm(const Vector& values) { + double norm = 0.0; + for (std::size_t i = 0U; i < values.Size(); ++i) + norm = (std::max)(norm, std::abs(values[i])); + return norm; +} +/// @brief Compares immutable semantic and raw source identities. +bool SameIdentity(const SourceEntityId& a, const SourceEntityId& b) { + return a.instance_name == b.instance_name && + a.source_label == b.source_label && + a.source_label_text == b.source_label_text; +} +/// @brief Builds one raw criterion with the approved branch and physical unit. +NewtonChannelDecision Decision(const char* name, double numerator, + double denominator, NewtonChannelUnit unit) { + const bool absolute = denominator == 0.0; + const double threshold = + absolute ? kNewtonZeroScaleThreshold : kNewtonRelativeThreshold; + return { + name, + numerator, + denominator, + absolute ? NewtonConvergenceBranch::kZeroScale + : NewtonConvergenceBranch::kRelative, + threshold, + absolute ? unit : NewtonChannelUnit::kDimensionless, + absolute ? numerator <= threshold : numerator / denominator <= threshold}; +} +} // namespace + +Result EvaluateSolidStaticConvergence( + const Vector& residual, const Vector& internal, const Vector& external, + const Vector& correction, const Vector& increment) { + const std::size_t size = residual.Size(); + if (!FiniteVector(residual, size) || !FiniteVector(internal, size) || + !FiniteVector(external, size) || !FiniteVector(correction, size) || + !FiniteVector(increment, size)) + return Result::Failure(StateFailure( + "nonfinite-solid-newton-data", + "Static channels require finite vectors of the same free dimension.")); + const double external_norm = InfinityNorm(external); + const double force_scale = + external_norm == 0.0 ? 0.0 + : (std::max)(external_norm, InfinityNorm(internal)); + SolidStaticIterationRecord record; + record.force_residual = Decision("residual_force", InfinityNorm(residual), + force_scale, NewtonChannelUnit::kForce); + record.translation_correction = + Decision("correction_translation", InfinityNorm(correction), + InfinityNorm(increment), NewtonChannelUnit::kLength); + record.convergence_evaluated = true; + return Result::Success(std::move(record)); +} + +Result SolidStaticState::Create(const Domain& domain, + const DofManager& dofs) { + if (!domain.IsSolid() || domain.Steps().Size() != 1U || + domain.Frequency() != nullptr || domain.Dynamic() != nullptr || + dofs.Profile() != DofProfile::kSolidStructural || + dofs.ComponentsPerNode() != 3U || + dofs.FullDofCount() != 3U * domain.Nodes().size()) + return Result::Failure( + StateFailure("invalid-solid-static-state", + "Static state requires one Solid static step and matching " + "3N mapping.")); + const Status mapping = dofs.ValidateInvariants(); + if (!mapping.IsOk()) return Result::Failure(mapping); + SolidStaticState state; + state.initialized_ = true; + state.nonlinear_ = domain.Steps()[0].IsNonlinearStatic(); + state.constrained_dofs_ = dofs.ConstrainedDofs(); + state.prescribed_values_ = dofs.PrescribedValues(); + const std::size_t count = dofs.FullDofCount(); + state.committed_.displacement = Vector{count}; + state.committed_.internal_force = Vector{count}; + state.committed_.external_force = Vector{count}; + state.committed_.residual = Vector{count}; + state.committed_.reaction = Vector{count}; + for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { + state.committed_.nodal_internal_order.push_back( + static_cast(node)); + state.committed_.nodal_source_ids.push_back(domain.Nodes()[node].source_id); + } + // Capture the immutable point coordinates from certified zero-state kernels. + // Later state commits cannot substitute another element or result location. + for (std::size_t e = 0U; e < domain.Elements().Size(); ++e) { + auto element = ElementFactory{}.Create(domain.Elements()[e], domain); + if (!element.HasValue()) + return Result::Failure(element.GetStatus()); + const auto* hex = dynamic_cast(element.Value().get()); + if (hex == nullptr) + return Result::Failure(StateFailure( + "invalid-solid-static-state", "Static state requires Hex8 kernels.")); + auto response = hex->EvaluateLinear(Vector{24U}); + if (!response.HasValue()) + return Result::Failure(response.GetStatus()); + state.reference_points_.insert(state.reference_points_.end(), + response.Value().points.begin(), + response.Value().points.end()); + } + return Result::Success(std::move(state)); +} + +Status SolidStaticState::CommitAcceptedFrame(const SolidStaticFrame& frame) { + const auto invalid = []() { + return StateFailure( + "invalid-solid-static-frame", + "Static candidate identity, projection or recovery is invalid."); + }; + if (!initialized_ || frame.frame_index != accepted_frames_.size() || + frame.increment_index == 0U || !std::isfinite(frame.load_factor) || + !(frame.load_factor > committed_.load_factor) || + frame.load_factor > 1.0 || + frame.is_factor_one != (frame.load_factor == 1.0) || + (nonlinear_ && + (frame.increment_index != committed_.increment_index + 1U || + frame.iteration_count == 0U || frame.iteration_count > 100U)) || + (!nonlinear_ && (!accepted_frames_.empty() || frame.load_factor != 1.0 || + frame.iteration_count != 0U))) + return invalid(); + const std::size_t count = committed_.displacement.Size(); + if (frame.nodal_internal_order != committed_.nodal_internal_order || + frame.nodal_source_ids.size() != committed_.nodal_source_ids.size() || + !FiniteVector(frame.displacement, count) || + !FiniteVector(frame.internal_force, count) || + !FiniteVector(frame.external_force, count) || + !FiniteVector(frame.residual, count) || + !FiniteVector(frame.reaction, count) || + frame.solid_results.size() != reference_points_.size() || + !std::isfinite(frame.physical_strain_energy) || + !Finite(frame.global_equilibrium) || + !Finite(frame.balance.internal_force_sum) || + !Finite(frame.balance.reference_internal_moment) || + !Finite(frame.balance.current_internal_moment) || + !Finite(frame.balance.free_residual_moment) || + !Finite(frame.balance.linearization_moment) || + !Finite(frame.balance.internal_force_scale) || + !Finite(frame.balance.reference_moment_scale) || + !Finite(frame.balance.current_moment_scale)) + return invalid(); + for (std::size_t node = 0U; node < frame.nodal_source_ids.size(); ++node) + if (!SameIdentity(frame.nodal_source_ids[node], + committed_.nodal_source_ids[node])) + return invalid(); + std::vector constrained(count, false); + for (std::size_t c = 0U; c < constrained_dofs_.size(); ++c) { + const auto full = constrained_dofs_[c]; + constrained[full] = true; + if (frame.displacement[full] != frame.load_factor * prescribed_values_[c]) + return invalid(); + } + for (std::size_t i = 0U; i < count; ++i) + if (frame.residual[i] != + frame.internal_force[i] - frame.external_force[i] || + frame.reaction[i] != (constrained[i] ? frame.residual[i] : 0.0)) + return invalid(); + for (std::size_t g = 0U; g < reference_points_.size(); ++g) { + const auto& row = frame.solid_results[g]; + const auto& reference = reference_points_[g]; + if (row.element != reference.element || + row.integration_point != reference.integration_point || + row.natural_coordinates != reference.natural_coordinates || + row.reference_coordinates != reference.reference_coordinates || + !Finite(row.cauchy_stress) || !Finite(row.strain)) + return invalid(); + } + // Prepare both copies before mutating committed history. + SolidStaticFrame next = frame; + accepted_frames_.push_back(frame); + committed_ = std::move(next); + return Status::Ok(); +} +const std::vector& SolidStaticState::AcceptedFrames() + const noexcept { + return accepted_frames_; +} +const SolidStaticFrame& SolidStaticState::CommittedState() const noexcept { + return committed_; +} +const std::vector& +SolidStaticState::IterationRecords() const noexcept { + return iterations_; +} +} // namespace fesa diff --git a/tests/unit/analysis/solid_static_procedure_test.cpp b/tests/unit/analysis/solid_static_procedure_test.cpp new file mode 100644 index 0000000..8567402 --- /dev/null +++ b/tests/unit/analysis/solid_static_procedure_test.cpp @@ -0,0 +1,738 @@ +#include "fesa/analysis/solid_static_procedure.h" + +#include + +#include +#include +#include +#include +#include +#include +#include + +#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 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::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 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::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::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(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::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, 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(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'; + } +} diff --git a/tests/unit/analysis/solid_static_state_test.cpp b/tests/unit/analysis/solid_static_state_test.cpp new file mode 100644 index 0000000..745713a --- /dev/null +++ b/tests/unit/analysis/solid_static_state_test.cpp @@ -0,0 +1,60 @@ +#include "fesa/analysis/solid_static_state.h" + +#include + +#include + +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::infinity(); + EXPECT_FALSE(fesa::EvaluateSolidStaticConvergence(r, f, p, du, u).HasValue()); + EXPECT_FALSE( + fesa::EvaluateSolidStaticConvergence(r, f, fesa::Vector{2U}, du, u) + .HasValue()); +}