From a5676814118bd56b506966bff40164a2c3355c3d Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Mon, 24 Aug 2026 14:13:14 +0900 Subject: [PATCH] feat(nonlinear-static-newton-raphson): step 5 - load-boundary-path-and-convergence-gates --- .../fesa/analysis/nonlinear_static_control.h | 131 +++++ src/fesa/CMakeLists.txt | 1 + .../analysis/nonlinear_static_control.cpp | 548 ++++++++++++++++++ tests/CMakeLists.txt | 1 + .../nonlinear_static_control_test.cpp | 365 ++++++++++++ 5 files changed, 1046 insertions(+) create mode 100644 include/fesa/analysis/nonlinear_static_control.h create mode 100644 src/fesa/analysis/nonlinear_static_control.cpp create mode 100644 tests/unit/analysis/nonlinear_static_control_test.cpp diff --git a/include/fesa/analysis/nonlinear_static_control.h b/include/fesa/analysis/nonlinear_static_control.h new file mode 100644 index 0000000..82a0be1 --- /dev/null +++ b/include/fesa/analysis/nonlinear_static_control.h @@ -0,0 +1,131 @@ +#ifndef FESA_ANALYSIS_NONLINEAR_STATIC_CONTROL_H_ +#define FESA_ANALYSIS_NONLINEAR_STATIC_CONTROL_H_ + +#include +#include +#include + +#include "fesa/core/diagnostic.h" +#include "fesa/core/status.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/math/sparse_matrix.h" +#include "fesa/math/vector.h" +#include "fesa/math/vector3.h" +#include "fesa/model/domain.h" + +namespace fesa { + +/// @brief Identifies one deterministic nonlinear load-factor increment. +struct LoadFactorIncrement { + std::uint32_t increment_index; + double load_factor; + bool is_factor_one; +}; + +/// @brief Identifies the physical dimension of a convergence threshold. +enum class NewtonChannelUnit { + kForce, + kForceLength, + kLength, + kDimensionless, +}; + +inline constexpr double kNewtonRelativeThreshold = 1.0e-2; +inline constexpr double kNewtonZeroScaleThreshold = 1.0e-3; + +/// @brief Carries one same-trial free residual convergence channel. +struct NewtonResidualChannelInput { + std::string identity; + Vector free_residual; + Vector free_internal; + Vector free_external; +}; + +/// @brief Carries one same-iteration free correction convergence channel. +struct NewtonCorrectionChannelInput { + std::string identity; + Vector fresh_correction; + Vector current_increment; +}; + +/// @brief Groups the four approved Newton convergence channels. +struct NewtonConvergenceInput { + NewtonResidualChannelInput force_residual; + NewtonResidualChannelInput moment_residual; + NewtonCorrectionChannelInput translation_correction; + NewtonCorrectionChannelInput rotation_correction; +}; + +/// @brief Identifies the only approved Newton convergence comparison branch. +enum class NewtonConvergenceBranch { kRelative, kZeroScale }; + +/// @brief Records the deterministic decision for one Newton channel. +struct NewtonChannelDecision { + std::string identity; + double numerator_norm; + double denominator_scale; + NewtonConvergenceBranch branch; + double threshold; + NewtonChannelUnit threshold_unit; + bool passed; +}; + +/// @brief Groups all four channel decisions and their logical conjunction. +struct NewtonConvergenceResult { + NewtonChannelDecision force_residual; + NewtonChannelDecision moment_residual; + NewtonChannelDecision translation_correction; + NewtonChannelDecision rotation_correction; + bool combined_passed; +}; + +/// @brief Stores the full-space fixed-load vector and its exact tangent. +struct GlobalFixedLoadContribution { + Vector generalized_external_force; + SparseMatrix generalized_external_tangent; + std::vector diagnostics; +}; + +/// @brief Derives the fixed nonlinear factor path from raw static controls. +/// @param step Immutable nonlinear static step that owns raw source controls. +/// @return Ordered factors ending in an exactly assigned factor one, or an +/// input failure. +Result> BuildFixedLoadFactorSequence( + const StepDefinition& step); + +/// @brief Scales factor-one prescribed targets from the zero step-start state. +/// @param step Immutable nonlinear step whose boundary order fixes output +/// order. +/// @param load_factor Current finite factor in the closed interval [0, 1]. +/// @return One absolute prescribed value per step boundary component, or an +/// input failure. +Result ApplyAbsolutePrescribedDisplacementTarget( + const StepDefinition& step, double load_factor); + +/// @brief Assembles fixed-global CLOAD work once in full generalized space. +/// @param domain Immutable nonlinear semantic model with one active step. +/// @param dofs Stable full-DOF mapping for domain nodes. +/// @param trial_generalized_coordinates Current full global [u, phi] trial. +/// @param load_factor Current finite factor in the closed interval [0, 1]. +/// @return External vector and tangent for r=f_int-f_ext and K=K_int-K_ext, +/// or a structured model/solver failure. +Result AssembleGlobalFixedLoadContribution( + const Domain& domain, const DofManager& dofs, + const Vector& trial_generalized_coordinates, double load_factor); + +/// @brief Rejects a nonzero source moment with director-parallel work. +/// @param director Current finite unit director in global coordinates. +/// @param source_moment Immutable fixed-global source moment. +/// @return Success when rho_M is at most 1e-12 or the source moment is zero. +Status CheckMomentAdmissibility(const Vector3& director, + const Vector3& source_moment); + +/// @brief Evaluates the four immutable, typed Newton convergence gates. +/// @param input Same-trial free residuals and same-iteration corrections. +/// @return Per-channel decision records and their logical conjunction. +NewtonConvergenceResult EvaluateNewtonConvergence( + const NewtonConvergenceInput& input); + +} // namespace fesa + +#endif // FESA_ANALYSIS_NONLINEAR_STATIC_CONTROL_H_ diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index de0913d..1f38875 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -4,6 +4,7 @@ add_library( analysis/analysis_model.cpp analysis/analysis_state.cpp analysis/linear_static_analysis.cpp + analysis/nonlinear_static_control.cpp app/fesa_application.cpp app/output_paths.cpp assembly/load_assembler.cpp diff --git a/src/fesa/analysis/nonlinear_static_control.cpp b/src/fesa/analysis/nonlinear_static_control.cpp new file mode 100644 index 0000000..543ce8f --- /dev/null +++ b/src/fesa/analysis/nonlinear_static_control.cpp @@ -0,0 +1,548 @@ +#include "fesa/analysis/nonlinear_static_control.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/assembly/load_assembler.h" +#include "fesa/elements/finite_rotation.h" + +namespace fesa { +namespace { + +constexpr std::size_t kDofsPerNode = 6U; +constexpr std::size_t kTranslationDofCount = 3U; +constexpr double kMomentAdmissibilityThreshold = 1.0e-12; + +/// @brief Creates a deterministic control-layer diagnostic status. +Status ControlFailure(const FailureCategory category, const std::string& code, + const SourceLocation& location, + const std::string& keyword, const std::string& identity, + const std::string& message) { + return Status::Failure(category, {{Severity::kError, code, location, keyword, + identity, message}}); +} + +/// @brief Verifies raw nonlinear static controls without deriving state. +Status ValidateNonlinearStaticControls(const StepDefinition& step) { + if (!step.IsNonlinearStatic() || !std::isfinite(step.InitialIncrement()) || + !(step.InitialIncrement() > 0.0) || !std::isfinite(step.TimePeriod()) || + step.TimePeriod() != 1.0 || !std::isfinite(step.MinimumIncrement()) || + !(step.MinimumIncrement() > 0.0) || + !std::isfinite(step.MaximumIncrement()) || + !(step.MaximumIncrement() > 0.0)) { + return ControlFailure( + FailureCategory::kInput, "invalid-nonlinear-static-data", + step.Location(), "STATIC", step.Name(), + "A nonlinear static step requires four positive finite controls and " + "an exactly factor-one time period."); + } + return Status::Ok(); +} + +/// @brief Returns a finite three-component vector from full trial storage. +Result RotationVectorForNode(const Vector& trial, + const std::size_t node) { + const std::size_t first = node * kDofsPerNode + kTranslationDofCount; + if (first + kTranslationDofCount > trial.Size()) { + return Result::Failure(ControlFailure( + FailureCategory::kSolver, "nonfinite-newton-data", {{}, 0U}, + "NONLINEAR_STATIC_CONTROL", std::to_string(node), + "The full trial coordinate vector does not contain every nodal " + "rotation component.")); + } + const Vector3 rotation{trial[first], trial[first + 1U], trial[first + 2U]}; + if (!rotation.IsFinite()) { + return Result::Failure(ControlFailure( + FailureCategory::kSolver, "nonfinite-newton-data", {{}, 0U}, + "NONLINEAR_STATIC_CONTROL", std::to_string(node), + "The current trial rotation coordinate must be finite.")); + } + return Result::Success(rotation); +} + +/// @brief Builds a proper initial frame from its semantic nodal record. +Matrix InitialFrame(const ShellNodeInitialFrame& frame) { + Matrix result{3U, 3U}; + for (std::size_t row = 0U; row < kTranslationDofCount; ++row) { + result(row, 0U) = frame.tangent_a[row]; + result(row, 1U) = frame.tangent_b[row]; + result(row, 2U) = frame.director[row]; + } + return result; +} + +/// @brief Multiplies a 3-by-3 matrix by a Vector3. +Vector3 Multiply(const Matrix& matrix, const Vector3& vector) { + return {matrix(0U, 0U) * vector.X() + matrix(0U, 1U) * vector.Y() + + matrix(0U, 2U) * vector.Z(), + matrix(1U, 0U) * vector.X() + matrix(1U, 1U) * vector.Y() + + matrix(1U, 2U) * vector.Z(), + matrix(2U, 0U) * vector.X() + matrix(2U, 1U) * vector.Y() + + matrix(2U, 2U) * vector.Z()}; +} + +/// @brief Multiplies a 3-by-3 matrix transpose by a Vector3. +Vector3 MultiplyTranspose(const Matrix& matrix, const Vector3& vector) { + return {matrix(0U, 0U) * vector.X() + matrix(1U, 0U) * vector.Y() + + matrix(2U, 0U) * vector.Z(), + matrix(0U, 1U) * vector.X() + matrix(1U, 1U) * vector.Y() + + matrix(2U, 1U) * vector.Z(), + matrix(0U, 2U) * vector.X() + matrix(1U, 2U) * vector.Y() + + matrix(2U, 2U) * vector.Z()}; +} + +/// @brief Forms the three-dimensional skew matrix of a global vector. +Matrix Skew(const Vector3& vector) { + Matrix result{3U, 3U}; + result(0U, 1U) = -vector.Z(); + result(0U, 2U) = vector.Y(); + result(1U, 0U) = vector.Z(); + result(1U, 2U) = -vector.X(); + result(2U, 0U) = -vector.Y(); + result(2U, 1U) = vector.X(); + return result; +} + +/// @brief Scales a small dense matrix without changing its shape. +Matrix Scale(const Matrix& matrix, const double scale) { + Matrix result{matrix.Rows(), matrix.Columns()}; + for (std::size_t row = 0U; row < matrix.Rows(); ++row) { + for (std::size_t column = 0U; column < matrix.Columns(); ++column) { + result(row, column) = scale * matrix(row, column); + } + } + return result; +} + +/// @brief Adds matching small dense matrices. +Matrix Add(const Matrix& left, const Matrix& right) { + Matrix result{left.Rows(), left.Columns()}; + for (std::size_t row = 0U; row < left.Rows(); ++row) { + for (std::size_t column = 0U; column < left.Columns(); ++column) { + result(row, column) = left(row, column) + right(row, column); + } + } + return result; +} + +/// @brief Evaluates left-Jacobian coefficients and derivatives by theta +/// squared. +void LeftJacobianCoefficients(const double theta_squared, double& a, double& b, + double& da_dtheta_squared, + double& db_dtheta_squared) { + if (theta_squared < 1.0e-8) { + const double theta_fourth = theta_squared * theta_squared; + const double theta_sixth = theta_fourth * theta_squared; + a = 0.5 - theta_squared / 24.0 + theta_fourth / 720.0 - + theta_sixth / 40320.0; + b = 1.0 / 6.0 - theta_squared / 120.0 + theta_fourth / 5040.0 - + theta_sixth / 362880.0; + da_dtheta_squared = + -1.0 / 24.0 + theta_squared / 360.0 - theta_fourth / 13440.0; + db_dtheta_squared = + -1.0 / 120.0 + theta_squared / 2520.0 - theta_fourth / 120960.0; + return; + } + + const double theta = std::sqrt(theta_squared); + const double sine = std::sin(theta); + const double cosine = std::cos(theta); + a = (1.0 - cosine) / theta_squared; + b = (theta - sine) / (theta_squared * theta); + da_dtheta_squared = (theta * sine - 2.0 * (1.0 - cosine)) / + (2.0 * theta_squared * theta_squared); + db_dtheta_squared = + (1.0 - cosine) / (2.0 * theta_squared * theta_squared) - + 3.0 * (theta - sine) / (2.0 * theta_squared * theta_squared * theta); +} + +/// @brief Differentiates the SO(3) left Jacobian along one coordinate vector. +Matrix LeftJacobianDerivative(const Vector3& phi, const Vector3& direction) { + const Matrix skew = Skew(phi); + const Matrix direction_skew = Skew(direction); + const Matrix skew_squared = skew.Multiply(skew); + double a = 0.0; + double b = 0.0; + double da_dtheta_squared = 0.0; + double db_dtheta_squared = 0.0; + LeftJacobianCoefficients(phi.Dot(phi), a, b, da_dtheta_squared, + db_dtheta_squared); + const double derivative_scale = 2.0 * phi.Dot(direction); + return Add(Add(Scale(skew, derivative_scale * da_dtheta_squared), + Scale(direction_skew, a)), + Add(Scale(skew_squared, derivative_scale * db_dtheta_squared), + Scale(Add(direction_skew.Multiply(skew), + skew.Multiply(direction_skew)), + b))); +} + +/// @brief Projects the immutable source moment into current director tangent +/// space. +Vector3 ProjectMoment(const Vector3& director, const Vector3& source_moment) { + return source_moment - director.Dot(source_moment) * director; +} + +/// @brief Differentiates the projected moment along a coordinate direction. +Vector3 ProjectMomentDerivative(const Vector3& director, + const Vector3& director_derivative, + const Vector3& source_moment) { + return -director.Dot(source_moment) * director_derivative - + director_derivative.Dot(source_moment) * director; +} + +/// @brief Computes the infinity norm and reports finite input validity. +double InfinityNorm(const Vector& values, bool& is_finite) { + double result = 0.0; + for (std::size_t index = 0U; index < values.Size(); ++index) { + if (!std::isfinite(values[index])) { + is_finite = false; + return std::numeric_limits::quiet_NaN(); + } + result = (std::max)(result, std::abs(values[index])); + } + return result; +} + +/// @brief Evaluates one residual channel with its fixed physical zero-scale +/// unit. +NewtonChannelDecision EvaluateResidualChannel( + const NewtonResidualChannelInput& input, const NewtonChannelUnit unit) { + bool is_finite = input.free_residual.Size() == input.free_internal.Size() && + input.free_residual.Size() == input.free_external.Size(); + const double numerator = InfinityNorm(input.free_residual, is_finite); + const double internal = InfinityNorm(input.free_internal, is_finite); + const double external = InfinityNorm(input.free_external, is_finite); + const double denominator = (std::max)(internal, external); + if (!is_finite) { + return {input.identity, + numerator, + denominator, + NewtonConvergenceBranch::kZeroScale, + kNewtonZeroScaleThreshold, + unit, + false}; + } + if (denominator == 0.0) { + return {input.identity, + numerator, + denominator, + NewtonConvergenceBranch::kZeroScale, + kNewtonZeroScaleThreshold, + unit, + numerator <= kNewtonZeroScaleThreshold}; + } + return {input.identity, + numerator, + denominator, + NewtonConvergenceBranch::kRelative, + kNewtonRelativeThreshold, + NewtonChannelUnit::kDimensionless, + numerator / denominator <= kNewtonRelativeThreshold}; +} + +/// @brief Evaluates one correction channel with its fixed physical zero-scale +/// unit. +NewtonChannelDecision EvaluateCorrectionChannel( + const NewtonCorrectionChannelInput& input, const NewtonChannelUnit unit) { + bool is_finite = + input.fresh_correction.Size() == input.current_increment.Size(); + const double numerator = InfinityNorm(input.fresh_correction, is_finite); + const double denominator = InfinityNorm(input.current_increment, is_finite); + if (!is_finite) { + return {input.identity, + numerator, + denominator, + NewtonConvergenceBranch::kZeroScale, + kNewtonZeroScaleThreshold, + unit, + false}; + } + if (denominator == 0.0) { + return {input.identity, + numerator, + denominator, + NewtonConvergenceBranch::kZeroScale, + kNewtonZeroScaleThreshold, + unit, + numerator <= kNewtonZeroScaleThreshold}; + } + return {input.identity, + numerator, + denominator, + NewtonConvergenceBranch::kRelative, + kNewtonRelativeThreshold, + NewtonChannelUnit::kDimensionless, + numerator / denominator <= kNewtonRelativeThreshold}; +} + +} // namespace + +Result> BuildFixedLoadFactorSequence( + const StepDefinition& step) { + const Status control_status = ValidateNonlinearStaticControls(step); + if (!control_status.IsOk()) { + return Result>::Failure(control_status); + } + + std::vector increments; + std::uint32_t increment_index = 1U; + while (true) { + const double ordinal_factor = + static_cast(increment_index) * step.InitialIncrement(); + if (!std::isfinite(ordinal_factor)) { + return Result>::Failure(ControlFailure( + FailureCategory::kInput, "invalid-nonlinear-static-data", + step.Location(), "STATIC", step.Name(), + "The fixed load-factor sequence must remain finite.")); + } + if (!(ordinal_factor < 1.0)) { + increments.push_back({increment_index, 1.0, true}); + break; + } + increments.push_back({increment_index, ordinal_factor, false}); + if (increment_index == (std::numeric_limits::max)()) { + return Result>::Failure(ControlFailure( + FailureCategory::kInput, "invalid-nonlinear-static-data", + step.Location(), "STATIC", step.Name(), + "The fixed load-factor sequence exceeds the supported increment " + "identity range.")); + } + ++increment_index; + } + return Result>::Success( + std::move(increments)); +} + +Result ApplyAbsolutePrescribedDisplacementTarget( + const StepDefinition& step, const double load_factor) { + const Status control_status = ValidateNonlinearStaticControls(step); + if (!control_status.IsOk()) { + return Result::Failure(control_status); + } + if (!std::isfinite(load_factor) || load_factor < 0.0 || load_factor > 1.0) { + return Result::Failure(ControlFailure( + FailureCategory::kInput, "invalid-nonlinear-static-data", + step.Location(), "BOUNDARY", step.Name(), + "The prescribed-displacement path requires a finite load factor in " + "the closed interval from zero to one.")); + } + + const auto& boundaries = step.PrescribedDisplacements(); + Vector targets{boundaries.Size()}; + for (std::size_t index = 0U; index < boundaries.Size(); ++index) { + const double target = boundaries[index].PrescribedValue(); + const double value = load_factor * target; + if (!std::isfinite(target) || !std::isfinite(value)) { + return Result::Failure(ControlFailure( + FailureCategory::kInput, "invalid-nonlinear-static-data", + boundaries[index].Location(), "BOUNDARY", + boundaries[index].Target().target_name_or_label, + "The factor-one prescribed displacement target must remain finite.")); + } + targets[index] = value; + } + return Result::Success(std::move(targets)); +} + +Status CheckMomentAdmissibility(const Vector3& director, + const Vector3& source_moment) { + if (!source_moment.IsFinite()) { + return ControlFailure(FailureCategory::kModel, "nonfinite-newton-data", + {{}, 0U}, "CLOAD", "", + "The fixed-global source moment must be finite."); + } + if (source_moment == Vector3{}) { + return Status::Ok(); + } + if (!director.IsFinite() || !(director.Norm() > 0.0)) { + return ControlFailure( + FailureCategory::kSolver, "nonfinite-newton-data", {{}, 0U}, "CLOAD", + "", "The current shell director must be finite and nonzero."); + } + const double source_norm = source_moment.Norm(); + const double rho = std::abs(director.Dot(source_moment)) / source_norm; + if (!std::isfinite(source_norm) || !(source_norm > 0.0) || + !std::isfinite(rho)) { + return ControlFailure( + FailureCategory::kModel, "nonfinite-newton-data", {{}, 0U}, "CLOAD", "", + "The fixed-global source moment norm must be finite."); + } + if (!(rho <= kMomentAdmissibilityThreshold)) { + return ControlFailure( + FailureCategory::kModel, "unsupported-drilling-load", {{}, 0U}, "CLOAD", + "", + "The aggregate nodal moment has an unsupported director-parallel " + "component."); + } + return Status::Ok(); +} + +Result AssembleGlobalFixedLoadContribution( + const Domain& domain, const DofManager& dofs, + const Vector& trial_generalized_coordinates, const double load_factor) { + if (!std::isfinite(load_factor) || load_factor < 0.0 || load_factor > 1.0 || + domain.Steps().Size() != 1U || !domain.Steps()[0U].IsNonlinearStatic()) { + return Result::Failure(ControlFailure( + FailureCategory::kModel, "invalid-nonlinear-static-data", + {domain.SourcePath(), 0U}, "NONLINEAR_STATIC_CONTROL", + domain.SourceContentIdentity(), + "Fixed-load assembly requires one nonlinear step and a finite load " + "factor in the closed interval from zero to one.")); + } + const Status dof_status = dofs.ValidateInvariants(); + if (!dof_status.IsOk()) { + return Result::Failure(dof_status); + } + if (dofs.FullDofCount() != trial_generalized_coordinates.Size() || + domain.Nodes().size() > + (std::numeric_limits::max)() / kDofsPerNode || + dofs.FullDofCount() != domain.Nodes().size() * kDofsPerNode) { + return Result::Failure(ControlFailure( + FailureCategory::kSolver, "nonfinite-newton-data", + {domain.SourcePath(), 0U}, "NONLINEAR_STATIC_CONTROL", + domain.SourceContentIdentity(), + "The full trial coordinates must match the stable Domain DOF " + "dimension.")); + } + for (std::size_t index = 0U; index < trial_generalized_coordinates.Size(); + ++index) { + if (!std::isfinite(trial_generalized_coordinates[index])) { + return Result::Failure(ControlFailure( + FailureCategory::kSolver, "nonfinite-newton-data", + {domain.SourcePath(), 0U}, "NONLINEAR_STATIC_CONTROL", + std::to_string(index), + "The full trial coordinate vector must contain finite values.")); + } + } + + auto model = AnalysisModel::Create(domain); + if (!model.HasValue()) { + return Result::Failure(model.GetStatus()); + } + auto source_load = LoadAssembler::AssembleFullNodalLoad(model.Value(), dofs); + if (!source_load.HasValue()) { + return Result::Failure( + source_load.GetStatus()); + } + + Vector external_force{dofs.FullDofCount()}; + for (std::size_t index = 0U; index < external_force.Size(); ++index) { + external_force[index] = load_factor * source_load.Value()[index]; + } + + std::vector frames(domain.Nodes().size(), + nullptr); + for (const auto& frame : domain.ShellNodeInitialFrames()) { + if (frame.node_index >= frames.size() || + frames[frame.node_index] != nullptr) { + return Result::Failure(ControlFailure( + FailureCategory::kModel, "invalid-shell-director", + {domain.SourcePath(), 0U}, "NODE", std::to_string(frame.node_index), + "Shell nodal frames must provide one unique record per node.")); + } + frames[frame.node_index] = &frame; + } + + std::vector tangent_contributions; + for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { + const std::size_t first_rotation = + node * kDofsPerNode + kTranslationDofCount; + const Vector3 source_moment{source_load.Value()[first_rotation], + source_load.Value()[first_rotation + 1U], + source_load.Value()[first_rotation + 2U]}; + if (source_moment == Vector3{}) { + continue; + } + if (frames[node] == nullptr) { + return Result::Failure(ControlFailure( + FailureCategory::kModel, "invalid-shell-director", + domain.Nodes()[node].location, "NODE", + domain.Nodes()[node].source_id.source_label_text, + "A node carrying a moment requires an initial shell frame.")); + } + auto phi = RotationVectorForNode(trial_generalized_coordinates, node); + if (!phi.HasValue()) { + return Result::Failure(phi.GetStatus()); + } + auto rotation = + MakeFiniteRotationState(phi.Value(), InitialFrame(*frames[node])); + if (!rotation.HasValue()) { + return Result::Failure(rotation.GetStatus()); + } + const Status admissibility = + CheckMomentAdmissibility(rotation.Value().director, source_moment); + if (!admissibility.IsOk()) { + return Result::Failure(admissibility); + } + + const Matrix left_jacobian = LeftJacobianSO3(phi.Value()); + const Vector3 projected = + ProjectMoment(rotation.Value().director, source_moment); + const Vector3 rotational_force = + load_factor * MultiplyTranspose(left_jacobian, projected); + for (std::size_t row = 0U; row < kTranslationDofCount; ++row) { + external_force[first_rotation + row] = rotational_force[row]; + } + + for (std::size_t column = 0U; column < kTranslationDofCount; ++column) { + Vector3 direction{}; + if (column == 0U) { + direction = {1.0, 0.0, 0.0}; + } else if (column == 1U) { + direction = {0.0, 1.0, 0.0}; + } else { + direction = {0.0, 0.0, 1.0}; + } + const Matrix left_jacobian_derivative = + LeftJacobianDerivative(phi.Value(), direction); + const Vector3 director_derivative = + Multiply(left_jacobian, direction).Cross(rotation.Value().director); + const Vector3 projected_derivative = ProjectMomentDerivative( + rotation.Value().director, director_derivative, source_moment); + const Vector3 tangent_column = + load_factor * + (MultiplyTranspose(left_jacobian_derivative, projected) + + MultiplyTranspose(left_jacobian, projected_derivative)); + for (std::size_t row = 0U; row < kTranslationDofCount; ++row) { + tangent_contributions.push_back( + {first_rotation + row, first_rotation + column, tangent_column[row], + node, row * kTranslationDofCount + column}); + } + } + } + + auto tangent = SparseMatrix::FromCoo(dofs.FullDofCount(), dofs.FullDofCount(), + std::move(tangent_contributions), + dofs.GetSparsePattern()); + if (!tangent.HasValue()) { + return Result::Failure(tangent.GetStatus()); + } + return Result::Success( + {std::move(external_force), std::move(tangent.Value()), {}}); +} + +NewtonConvergenceResult EvaluateNewtonConvergence( + const NewtonConvergenceInput& input) { + NewtonConvergenceResult result{ + EvaluateResidualChannel(input.force_residual, NewtonChannelUnit::kForce), + EvaluateResidualChannel(input.moment_residual, + NewtonChannelUnit::kForceLength), + EvaluateCorrectionChannel(input.translation_correction, + NewtonChannelUnit::kLength), + EvaluateCorrectionChannel(input.rotation_correction, + NewtonChannelUnit::kDimensionless), + false}; + result.combined_passed = + result.force_residual.passed && result.moment_residual.passed && + result.translation_correction.passed && result.rotation_correction.passed; + return result; +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d94a2a1..e39e5be 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -6,6 +6,7 @@ add_executable( unit/analysis/analysis_test.cpp unit/analysis/analysis_model_test.cpp unit/analysis/analysis_state_test.cpp + unit/analysis/nonlinear_static_control_test.cpp unit/assembly/parallel_for_test.cpp unit/assembly/load_assembler_test.cpp unit/assembly/sparse_assembler_test.cpp diff --git a/tests/unit/analysis/nonlinear_static_control_test.cpp b/tests/unit/analysis/nonlinear_static_control_test.cpp new file mode 100644 index 0000000..a64c959 --- /dev/null +++ b/tests/unit/analysis/nonlinear_static_control_test.cpp @@ -0,0 +1,365 @@ +#include "fesa/analysis/nonlinear_static_control.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/model/domain.h" +#include "fesa/model/model_types.h" + +namespace { + +fesa::Vector VectorFrom(const std::vector& values) { + fesa::Vector result{values.size()}; + for (std::size_t index = 0U; index < values.size(); ++index) { + result[index] = values[index]; + } + return result; +} + +struct Fixture { + std::unique_ptr domain; + std::unique_ptr model; + std::unique_ptr dofs; +}; + +Fixture MakeFixture( + std::vector boundaries = {}, + std::vector loads = {}, + const double initial_increment = 0.1, const double time_period = 1.0, + const double minimum_increment = 1.0e-9, + const double maximum_increment = 5.0) { + const std::filesystem::path source{"nonlinear-static-control.inp"}; + fesa::ModelDefinition definition{}; + definition.source_path = source; + definition.source_content_identity = "control-fixture"; + definition.nodes = { + {{"Shell-1", 10, "10"}, {0.0, 0.0, 0.0}, {source, 1U}}, + {{"Shell-1", 20, "20"}, {1.0, 0.0, 0.0}, {source, 2U}}, + {{"Shell-1", 30, "30"}, {1.0, 1.0, 0.0}, {source, 3U}}, + {{"Shell-1", 40, "40"}, {0.0, 1.0, 0.0}, {source, 4U}}, + {{"Shell-1", 50, "50"}, {-1.0, 1.0, 0.0}, {source, 5U}}, + {{"Shell-1", 60, "60"}, {-1.0, 0.0, 0.0}, {source, 6U}}, + }; + definition.materials = {{"Material", 1000.0, 0.25, {source, 7U}}}; + definition.shell_sections = {{"Shell", 0.1, 0U, {source, 8U}}}; + definition.shell_elements = { + {{"Shell-1", 1, "1"}, + fesa::ShellSourceElementType::kS4, + {0U, 1U, 2U, 3U}, + 0U, + 0U, + {source, 9U}}, + {{"Shell-1", 2, "2"}, + fesa::ShellSourceElementType::kS4, + {0U, 3U, 4U, 5U}, + 0U, + 0U, + {source, 10U}}, + }; + for (std::size_t node = 0U; node < definition.nodes.size(); ++node) { + definition.shell_node_initial_frames.push_back( + {static_cast(node), + {0.0, 0.0, 1.0}, + {1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}}); + } + definition.steps = {{"Step-1", + std::move(boundaries), + std::move(loads), + fesa::StaticProcedureKind::kNonlinearFullNewton, + true, + initial_increment, + time_period, + minimum_increment, + maximum_increment, + {source, 11U}}}; + + auto domain_result = fesa::Domain::Create(std::move(definition)); + if (!domain_result.HasValue()) { + throw std::runtime_error{"Control fixture Domain construction failed."}; + } + auto domain = + std::make_unique(std::move(domain_result.Value())); + auto model_result = fesa::AnalysisModel::Create(*domain); + if (!model_result.HasValue()) { + throw std::runtime_error{ + "Control fixture AnalysisModel construction failed."}; + } + auto model = + std::make_unique(std::move(model_result.Value())); + auto dofs_result = fesa::DofManager::Create(*model); + if (!dofs_result.HasValue()) { + throw std::runtime_error{"Control fixture DofManager construction failed."}; + } + auto dofs = + std::make_unique(std::move(dofs_result.Value())); + return {std::move(domain), std::move(model), std::move(dofs)}; +} + +fesa::Vector TrialCoordinates(const Fixture& fixture) { + return fesa::Vector{fixture.dofs->FullDofCount()}; +} + +fesa::Vector AddScaled(const fesa::Vector& values, + const fesa::Vector& direction, const double scale) { + fesa::Vector result = values; + result.Axpy(scale, direction); + return result; +} + +double MaximumAbsoluteDifference(const fesa::Vector& left, + const fesa::Vector& right) { + if (left.Size() != right.Size()) { + throw std::invalid_argument{"Vector dimensions must match."}; + } + double maximum = 0.0; + for (std::size_t index = 0U; index < left.Size(); ++index) { + maximum = (std::max)(maximum, std::abs(left[index] - right[index])); + } + return maximum; +} + +fesa::NewtonConvergenceInput RelativeConvergenceInput() { + return { + {"residual_force", VectorFrom({1.0e-2}), VectorFrom({1.0}), + VectorFrom({0.5})}, + {"residual_moment", VectorFrom({1.0e-2}), VectorFrom({1.0}), + VectorFrom({0.25})}, + {"correction_translation", VectorFrom({1.0e-2}), VectorFrom({1.0})}, + {"correction_rotation", VectorFrom({1.0e-2}), VectorFrom({1.0})}, + }; +} + +// NSNR-CTRL-001, NSNR-CTRL-002, and NSNR-CTRL-017 +TEST(NonlinearStaticControl, BuildsFixedFactorSequenceWithExactEndpoint) { + const auto fixture = MakeFixture(); + const auto result = + fesa::BuildFixedLoadFactorSequence(fixture.domain->Steps()[0U]); + + ASSERT_TRUE(result.HasValue()); + ASSERT_EQ(result.Value().size(), 10U); + for (std::size_t index = 0U; index < result.Value().size(); ++index) { + const auto& increment = result.Value()[index]; + EXPECT_EQ(increment.increment_index, + static_cast(index + 1U)); + EXPECT_DOUBLE_EQ(increment.load_factor, + static_cast(index + 1U) * 0.1); + EXPECT_EQ(increment.is_factor_one, index + 1U == result.Value().size()); + } + EXPECT_DOUBLE_EQ(result.Value().back().load_factor, 1.0); + + const auto one_hundred = MakeFixture({}, {}, 0.01); + const auto one_hundred_result = + fesa::BuildFixedLoadFactorSequence(one_hundred.domain->Steps()[0U]); + ASSERT_TRUE(one_hundred_result.HasValue()); + ASSERT_EQ(one_hundred_result.Value().size(), 100U); + EXPECT_EQ(one_hundred_result.Value().back().increment_index, 100U); + EXPECT_TRUE(one_hundred_result.Value().back().is_factor_one); +} + +// NSNR-CTRL-003 +TEST(NonlinearStaticControl, IgnoresPositiveMinimumAndMaximumDuringStepping) { + const auto first = MakeFixture({}, {}, 0.3, 1.0, 1.0e-12, 100.0); + const auto second = MakeFixture({}, {}, 0.3, 1.0, 0.29, 0.31); + const auto first_result = + fesa::BuildFixedLoadFactorSequence(first.domain->Steps()[0U]); + const auto second_result = + fesa::BuildFixedLoadFactorSequence(second.domain->Steps()[0U]); + + ASSERT_TRUE(first_result.HasValue()); + ASSERT_TRUE(second_result.HasValue()); + ASSERT_EQ(first_result.Value().size(), second_result.Value().size()); + for (std::size_t index = 0U; index < first_result.Value().size(); ++index) { + EXPECT_DOUBLE_EQ(first_result.Value()[index].load_factor, + second_result.Value()[index].load_factor); + EXPECT_EQ(first_result.Value()[index].increment_index, + second_result.Value()[index].increment_index); + } +} + +// NSNR-CTRL-004 and NSNR-CTRL-005 +TEST(NonlinearLoadBoundary, AppliesAbsoluteBoundaryTargetsWithoutCorrection) { + const std::filesystem::path source{"nonlinear-static-control.inp"}; + const auto fixture = MakeFixture( + {{"10", 1, 1, 4.0, {source, 12U}}, {"20", 4, 4, -0.5, {source, 13U}}}); + const auto result = fesa::ApplyAbsolutePrescribedDisplacementTarget( + fixture.domain->Steps()[0U], 0.25); + + ASSERT_TRUE(result.HasValue()); + ASSERT_EQ(result.Value().Size(), 2U); + EXPECT_DOUBLE_EQ(result.Value()[0U], 1.0); + EXPECT_DOUBLE_EQ(result.Value()[1U], -0.125); + for (const std::size_t constrained : fixture.dofs->ConstrainedDofs()) { + EXPECT_EQ(TrialCoordinates(fixture)[constrained], 0.0); + } +} + +// NSNR-CTRL-006 and NSNR-CTRL-008 +TEST(NonlinearLoadBoundary, AssemblesEverySharedNodeLoadOnceInFixedGlobalAxes) { + const std::filesystem::path source{"nonlinear-static-control.inp"}; + const std::vector loads{ + {"10", 1, 8.0, {source, 20U}}, + {"10", 2, -4.0, {source, 21U}}, + {"10", 4, 2.0, {source, 22U}}, + }; + const auto fixture = MakeFixture({}, loads); + const fesa::Vector trial = TrialCoordinates(fixture); + const auto result = fesa::AssembleGlobalFixedLoadContribution( + *fixture.domain, *fixture.dofs, trial, 0.25); + + ASSERT_TRUE(result.HasValue()); + ASSERT_EQ(result.Value().generalized_external_force.Size(), + fixture.dofs->FullDofCount()); + EXPECT_DOUBLE_EQ(result.Value().generalized_external_force[0U], 2.0); + EXPECT_DOUBLE_EQ(result.Value().generalized_external_force[1U], -1.0); + EXPECT_DOUBLE_EQ(result.Value().generalized_external_force[3U], 0.5); + EXPECT_DOUBLE_EQ(fixture.domain->Steps()[0U] + .ConcentratedLoads()[0U] + .GlobalComponents()[0U], + 8.0); + EXPECT_DOUBLE_EQ(fixture.domain->Steps()[0U] + .ConcentratedLoads()[2U] + .GlobalComponents()[3U], + 2.0); +} + +// NSNR-CTRL-007 +TEST(NonlinearLoadBoundary, RejectsDirectorAlignedMomentAboveTheExactBand) { + const fesa::Vector3 director{0.0, 0.0, 1.0}; + EXPECT_TRUE( + fesa::CheckMomentAdmissibility(director, {1.0, 0.0, 1.0e-12}).IsOk()); + const auto rejected = + fesa::CheckMomentAdmissibility(director, {1.0, 0.0, 2.0e-12}); + EXPECT_FALSE(rejected.IsOk()); + ASSERT_EQ(rejected.Diagnostics().size(), 1U); + EXPECT_EQ(rejected.Diagnostics()[0U].code, "unsupported-drilling-load"); +} + +// NSNR-CTRL-009 +TEST(NonlinearLoadBoundary, PreservesSourceOrderDuringGlobalAccumulation) { + const std::filesystem::path source{"nonlinear-static-control.inp"}; + const auto fixture = MakeFixture({}, {{"10", 1, 1.0e16, {source, 30U}}, + {"10", 1, -1.0e16, {source, 31U}}, + {"10", 1, 1.0, {source, 32U}}}); + const fesa::Vector trial = TrialCoordinates(fixture); + const auto first = fesa::AssembleGlobalFixedLoadContribution( + *fixture.domain, *fixture.dofs, trial, 1.0); + const auto second = fesa::AssembleGlobalFixedLoadContribution( + *fixture.domain, *fixture.dofs, trial, 1.0); + + ASSERT_TRUE(first.HasValue()); + ASSERT_TRUE(second.HasValue()); + EXPECT_DOUBLE_EQ(first.Value().generalized_external_force[0U], 1.0); + EXPECT_EQ(first.Value().generalized_external_tangent.RowOffsets(), + second.Value().generalized_external_tangent.RowOffsets()); + EXPECT_EQ(first.Value().generalized_external_tangent.ColumnIndices(), + second.Value().generalized_external_tangent.ColumnIndices()); + EXPECT_EQ(first.Value().generalized_external_tangent.Values(), + second.Value().generalized_external_tangent.Values()); +} + +// NSNR-CTRL-010 and NSNR-CTRL-011 +TEST(NonlinearLoadBoundary, + DifferentiatesOneFixedSpatialMomentAndUsesSubtraction) { + const std::filesystem::path source{"nonlinear-static-control.inp"}; + const auto fixture = MakeFixture({}, {{"10", 4, 3.0, {source, 40U}}}); + fesa::Vector trial = TrialCoordinates(fixture); + trial[5U] = 0.4; + const auto result = fesa::AssembleGlobalFixedLoadContribution( + *fixture.domain, *fixture.dofs, trial, 0.7); + ASSERT_TRUE(result.HasValue()); + + fesa::Vector direction{fixture.dofs->FullDofCount()}; + direction[5U] = 1.0; + constexpr double kStep = 1.0e-6; + const auto plus = fesa::AssembleGlobalFixedLoadContribution( + *fixture.domain, *fixture.dofs, AddScaled(trial, direction, kStep), 0.7); + const auto minus = fesa::AssembleGlobalFixedLoadContribution( + *fixture.domain, *fixture.dofs, AddScaled(trial, direction, -kStep), 0.7); + ASSERT_TRUE(plus.HasValue()); + ASSERT_TRUE(minus.HasValue()); + + fesa::Vector derivative = plus.Value().generalized_external_force; + derivative.Axpy(-1.0, minus.Value().generalized_external_force); + derivative.Scale(0.5 / kStep); + const fesa::Vector tangent_action = + result.Value().generalized_external_tangent.Multiply(direction); + EXPECT_LE(MaximumAbsoluteDifference(derivative, tangent_action), 1.0e-8); + + fesa::Vector internal{fixture.dofs->FullDofCount()}; + internal[3U] = 5.0; + fesa::Vector residual = internal; + residual.Axpy(-1.0, result.Value().generalized_external_force); + EXPECT_DOUBLE_EQ(residual[3U], + 5.0 - result.Value().generalized_external_force[3U]); + EXPECT_DOUBLE_EQ( + -tangent_action[3U], + -result.Value().generalized_external_tangent.Multiply(direction)[3U]); +} + +// NSNR-CTRL-012, NSNR-CTRL-013, and NSNR-CTRL-014 +TEST(NonlinearConvergence, UsesOwnedRelativeThresholdsForEveryChannel) { + static_assert(fesa::kNewtonRelativeThreshold == 1.0e-2); + static_assert(fesa::kNewtonZeroScaleThreshold == 1.0e-3); + const auto result = + fesa::EvaluateNewtonConvergence(RelativeConvergenceInput()); + + EXPECT_TRUE(result.combined_passed); + const std::array decisions{ + result.force_residual, result.moment_residual, + result.translation_correction, result.rotation_correction}; + for (const auto& decision : decisions) { + EXPECT_FALSE(decision.identity.empty()); + EXPECT_DOUBLE_EQ(decision.numerator_norm, 1.0e-2); + EXPECT_DOUBLE_EQ(decision.denominator_scale, 1.0); + EXPECT_EQ(decision.branch, fesa::NewtonConvergenceBranch::kRelative); + EXPECT_DOUBLE_EQ(decision.threshold, fesa::kNewtonRelativeThreshold); + EXPECT_EQ(decision.threshold_unit, fesa::NewtonChannelUnit::kDimensionless); + EXPECT_TRUE(decision.passed); + } +} + +// NSNR-CTRL-015 and NSNR-CTRL-016 +TEST(NonlinearConvergence, AppliesTypedZeroScaleGatesToEmptyChannels) { + const fesa::NewtonConvergenceInput input{ + {"residual_force", fesa::Vector{0U}, fesa::Vector{0U}, fesa::Vector{0U}}, + {"residual_moment", fesa::Vector{0U}, fesa::Vector{0U}, fesa::Vector{0U}}, + {"correction_translation", fesa::Vector{0U}, fesa::Vector{0U}}, + {"correction_rotation", fesa::Vector{0U}, fesa::Vector{0U}}, + }; + const auto result = fesa::EvaluateNewtonConvergence(input); + const std::array decisions{ + result.force_residual, result.moment_residual, + result.translation_correction, result.rotation_correction}; + const std::array units{ + fesa::NewtonChannelUnit::kForce, fesa::NewtonChannelUnit::kForceLength, + fesa::NewtonChannelUnit::kLength, + fesa::NewtonChannelUnit::kDimensionless}; + + EXPECT_TRUE(result.combined_passed); + for (std::size_t index = 0U; index < decisions.size(); ++index) { + EXPECT_DOUBLE_EQ(decisions[index].numerator_norm, 0.0); + EXPECT_DOUBLE_EQ(decisions[index].denominator_scale, 0.0); + EXPECT_EQ(decisions[index].branch, + fesa::NewtonConvergenceBranch::kZeroScale); + EXPECT_DOUBLE_EQ(decisions[index].threshold, + fesa::kNewtonZeroScaleThreshold); + EXPECT_EQ(decisions[index].threshold_unit, units[index]); + EXPECT_TRUE(decisions[index].passed); + } +} + +} // namespace