From f0f38b98039727948655e9eb149b2eb9e6d38d52 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Mon, 24 Aug 2026 15:49:08 +0900 Subject: [PATCH] feat(nonlinear-static-newton-raphson): step 10 - prepared-analysis-flow-and-full-procedure --- .../fesa/analysis/linear_static_analysis.h | 46 +- .../fesa/analysis/linear_static_procedure.h | 31 + .../analysis/nonlinear_static_procedure.h | 33 + include/fesa/analysis/prepared_analysis.h | 71 ++ src/fesa/CMakeLists.txt | 3 + src/fesa/analysis/linear_static_analysis.cpp | 183 +--- src/fesa/analysis/linear_static_procedure.cpp | 105 +++ .../analysis/nonlinear_static_procedure.cpp | 807 ++++++++++++++++++ src/fesa/analysis/prepared_analysis.cpp | 79 ++ src/fesa/app/fesa_application.cpp | 63 +- tests/CMakeLists.txt | 1 + .../analysis/linear_static_analysis_test.cpp | 1 + .../nonlinear_static_procedure_test.cpp | 157 ++++ .../integration/app/fesa_application_test.cpp | 14 + .../b33_reference_comparison_test.cpp | 1 + .../reference/mitc4_reference_cases_test.cpp | 1 + 16 files changed, 1373 insertions(+), 223 deletions(-) create mode 100644 include/fesa/analysis/linear_static_procedure.h create mode 100644 include/fesa/analysis/nonlinear_static_procedure.h create mode 100644 include/fesa/analysis/prepared_analysis.h create mode 100644 src/fesa/analysis/linear_static_procedure.cpp create mode 100644 src/fesa/analysis/nonlinear_static_procedure.cpp create mode 100644 src/fesa/analysis/prepared_analysis.cpp create mode 100644 tests/integration/analysis/nonlinear_static_procedure_test.cpp diff --git a/include/fesa/analysis/linear_static_analysis.h b/include/fesa/analysis/linear_static_analysis.h index 8911bde..f0e7595 100644 --- a/include/fesa/analysis/linear_static_analysis.h +++ b/include/fesa/analysis/linear_static_analysis.h @@ -1,18 +1,7 @@ #ifndef FESA_ANALYSIS_LINEAR_STATIC_ANALYSIS_H_ #define FESA_ANALYSIS_LINEAR_STATIC_ANALYSIS_H_ -#include -#include - #include "fesa/analysis/analysis.h" -#include "fesa/analysis/analysis_model.h" -#include "fesa/analysis/analysis_state.h" -#include "fesa/constraints/essential_constraint_policy.h" -#include "fesa/elements/element.h" -#include "fesa/fem/dof_manager.h" -#include "fesa/math/sparse_matrix.h" -#include "fesa/math/vector.h" -#include "fesa/model/domain.h" namespace fesa { @@ -20,9 +9,8 @@ class LinearSolver; class ParallelFor; class ResultsWriter; -/// @brief Orchestrates the single-step linear-static procedure. -/// @note Factorization, substitution, recovery, and writing remain separately -/// observable through injected backend boundaries. +/// @brief Retains the legacy Analysis facade for linear-static callers. +/// @note Numerical stages are owned by LinearStaticProcedureExecutor. class LinearStaticAnalysis final : public Analysis { public: /// @brief Creates a procedure using non-owned backend adapters. @@ -31,42 +19,14 @@ class LinearStaticAnalysis final : public Analysis { LinearSolver& linear_solver, ResultsWriter& results_writer); - /// @brief Runs the approved eight-stage linear-static lifecycle. + /// @brief Prepares, executes, and writes a linear-static analysis. /// @return The first stage failure or successful result finalization. Status Run(const AnalysisRequest& request) override; private: - /// @brief Initializes owned input and Domain state for a run candidate. - Status InitializeCandidate(const AnalysisRequest& request); - /// @brief Builds the non-owning active-model view. - Status BuildAnalysisModel(); - /// @brief Creates runtime elements, stable DOFs, and the sparse pattern. - Status BuildDofMapAndSparsePattern(); - /// @brief Assembles full stiffness and stable constraint partitions. - Status AssembleAndPartitionStiffness(); - /// @brief Factorizes Kff before any load assembly. - Status FactorizeFreeSystem(); - /// @brief Assembles loads and forms Ff-Kfc*dc without solving. - Status AssembleLoadsAndEffectiveRhs(); - /// @brief Substitutes the retained factorization and reconstructs full d. - Status SubstituteAndReconstruct(); - /// @brief Recovers a complete candidate before writing authoritative output. - Status RecoverAndWrite(); - const ParallelFor& parallel_for_; LinearSolver& linear_solver_; ResultsWriter& results_writer_; - AnalysisRequest request_; - std::unique_ptr domain_; - std::unique_ptr model_; - std::vector> elements_; - ElementView element_view_; - std::unique_ptr dofs_; - std::unique_ptr state_; - std::unique_ptr full_stiffness_; - std::unique_ptr partitioned_stiffness_; - std::unique_ptr effective_rhs_; - std::vector diagnostics_; }; } // namespace fesa diff --git a/include/fesa/analysis/linear_static_procedure.h b/include/fesa/analysis/linear_static_procedure.h new file mode 100644 index 0000000..fe02c04 --- /dev/null +++ b/include/fesa/analysis/linear_static_procedure.h @@ -0,0 +1,31 @@ +#ifndef FESA_ANALYSIS_LINEAR_STATIC_PROCEDURE_H_ +#define FESA_ANALYSIS_LINEAR_STATIC_PROCEDURE_H_ + +#include "fesa/analysis/prepared_analysis.h" + +namespace fesa { + +class LinearSolver; +class ParallelFor; + +/// @brief Executes the approved linear-static procedure on a prepared Domain. +class LinearStaticProcedureExecutor { + public: + /// @brief Creates a procedure executor with caller-owned backend adapters. + /// @param parallel_for Deterministic element-local execution adapter. + /// @param linear_solver Factorization and substitution adapter. + LinearStaticProcedureExecutor(const ParallelFor& parallel_for, + LinearSolver& linear_solver); + + /// @brief Computes linear-static state without writing HDF5 output. + /// @return A validated state and sorted diagnostics, or a procedure failure. + Result Execute(const Domain& domain) const; + + private: + const ParallelFor& parallel_for_; + LinearSolver& linear_solver_; +}; + +} // namespace fesa + +#endif // FESA_ANALYSIS_LINEAR_STATIC_PROCEDURE_H_ diff --git a/include/fesa/analysis/nonlinear_static_procedure.h b/include/fesa/analysis/nonlinear_static_procedure.h new file mode 100644 index 0000000..ba12a38 --- /dev/null +++ b/include/fesa/analysis/nonlinear_static_procedure.h @@ -0,0 +1,33 @@ +#ifndef FESA_ANALYSIS_NONLINEAR_STATIC_PROCEDURE_H_ +#define FESA_ANALYSIS_NONLINEAR_STATIC_PROCEDURE_H_ + +#include "fesa/analysis/prepared_analysis.h" + +namespace fesa { + +class NewtonIterationLogWriter; +class ParallelFor; +class TangentSolver; + +/// @brief Executes fixed-load-factor full-Newton nonlinear static analysis. +class NonlinearStaticProcedureExecutor { + public: + /// @brief Creates a procedure executor with caller-owned backend adapters. + /// @param parallel_for Deterministic element-local execution adapter. + /// @param tangent_solver Fresh unsymmetric tangent solve adapter. + NonlinearStaticProcedureExecutor(const ParallelFor& parallel_for, + TangentSolver& tangent_solver); + + /// @brief Computes accepted nonlinear frames and appends iteration records. + /// @return A validated nonlinear state without writing HDF5 output. + Result Execute( + const Domain& domain, NewtonIterationLogWriter& log) const; + + private: + const ParallelFor& parallel_for_; + TangentSolver& tangent_solver_; +}; + +} // namespace fesa + +#endif // FESA_ANALYSIS_NONLINEAR_STATIC_PROCEDURE_H_ diff --git a/include/fesa/analysis/prepared_analysis.h b/include/fesa/analysis/prepared_analysis.h new file mode 100644 index 0000000..22d7b4b --- /dev/null +++ b/include/fesa/analysis/prepared_analysis.h @@ -0,0 +1,71 @@ +#ifndef FESA_ANALYSIS_PREPARED_ANALYSIS_H_ +#define FESA_ANALYSIS_PREPARED_ANALYSIS_H_ + +#include +#include +#include + +#include "fesa/analysis/analysis_state.h" +#include "fesa/analysis/nonlinear_analysis_state.h" +#include "fesa/core/diagnostic.h" +#include "fesa/core/status.h" +#include "fesa/model/domain.h" + +namespace fesa { + +class LinearStaticProcedureExecutor; +class NewtonIterationLogWriter; +class NonlinearResultsWriter; +class NonlinearStaticProcedureExecutor; +class ResultsWriter; +struct DerivedOutputPaths; + +/// @brief Owns one parsed and semantically mapped analysis input. +struct PreparedAnalysisInput { + Domain domain; + std::vector diagnostics; +}; + +/// @brief Parses and maps an input file exactly once for procedure execution. +/// @param input_path Source Abaqus input path. +/// @return A movable Domain and sorted mapping diagnostics, or an input +/// failure. +Result PrepareAnalysisInput( + const std::filesystem::path& input_path); + +/// @brief Stores a completed linear-static calculation before result writing. +struct LinearStaticExecutionResult { + AnalysisState state; + std::vector diagnostics; +}; + +/// @brief Stores a completed nonlinear-static calculation before result +/// writing. +struct NonlinearStaticExecutionResult { + NonlinearAnalysisState state; + std::vector diagnostics; +}; + +using PreparedAnalysisOutcome = + std::variant; + +/// @brief Selects the approved procedure from a prepared Domain. +/// @pre The executors and log outlive this call. +/// @return A validated procedure outcome without opening a result file. +Result ExecutePreparedAnalysis( + const Domain& domain, LinearStaticProcedureExecutor& linear_executor, + NonlinearStaticProcedureExecutor& nonlinear_executor, + NewtonIterationLogWriter& log); + +/// @brief Writes an already computed outcome to its authoritative result file. +/// @pre The Newton log has closed successfully before this function is called. +/// @return Success after the selected writer atomically finalizes its result. +Status FinalizePreparedAnalysisResult( + const DerivedOutputPaths& paths, const Domain& domain, + const PreparedAnalysisOutcome& outcome, + const std::vector& diagnostics, ResultsWriter& linear_writer, + NonlinearResultsWriter& nonlinear_writer); + +} // namespace fesa + +#endif // FESA_ANALYSIS_PREPARED_ANALYSIS_H_ diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index 7a1afe6..5de593c 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -4,9 +4,12 @@ add_library( analysis/analysis_model.cpp analysis/analysis_state.cpp analysis/linear_static_analysis.cpp + analysis/linear_static_procedure.cpp analysis/newton_iteration_log.cpp analysis/nonlinear_analysis_state.cpp analysis/nonlinear_static_control.cpp + analysis/nonlinear_static_procedure.cpp + analysis/prepared_analysis.cpp app/fesa_application.cpp app/output_paths.cpp assembly/load_assembler.cpp diff --git a/src/fesa/analysis/linear_static_analysis.cpp b/src/fesa/analysis/linear_static_analysis.cpp index 80a6524..8377207 100644 --- a/src/fesa/analysis/linear_static_analysis.cpp +++ b/src/fesa/analysis/linear_static_analysis.cpp @@ -2,50 +2,12 @@ #include -#include "fesa/assembly/load_assembler.h" -#include "fesa/assembly/parallel_for.h" -#include "fesa/assembly/sparse_assembler.h" -#include "fesa/elements/element_factory.h" -#include "fesa/io/abaqus/domain_mapper.h" -#include "fesa/io/abaqus/input_reader.h" -#include "fesa/results/result_recovery.h" +#include "fesa/analysis/linear_static_procedure.h" +#include "fesa/analysis/prepared_analysis.h" #include "fesa/results/results_writer.h" -#include "fesa/solvers/linear/linear_solver.h" namespace fesa { -Status LinearStaticAnalysis::Run(const AnalysisRequest& request) { - Status status = InitializeCandidate(request); - if (!status.IsOk()) { - return status; - } - status = BuildAnalysisModel(); - if (!status.IsOk()) { - return status; - } - status = BuildDofMapAndSparsePattern(); - if (!status.IsOk()) { - return status; - } - status = AssembleAndPartitionStiffness(); - if (!status.IsOk()) { - return status; - } - status = FactorizeFreeSystem(); - if (!status.IsOk()) { - return status; - } - status = AssembleLoadsAndEffectiveRhs(); - if (!status.IsOk()) { - return status; - } - status = SubstituteAndReconstruct(); - if (!status.IsOk()) { - return status; - } - return RecoverAndWrite(); -} - LinearStaticAnalysis::LinearStaticAnalysis(const ParallelFor& parallel_for, LinearSolver& linear_solver, ResultsWriter& results_writer) @@ -53,138 +15,23 @@ LinearStaticAnalysis::LinearStaticAnalysis(const ParallelFor& parallel_for, linear_solver_{linear_solver}, results_writer_{results_writer} {} -Status LinearStaticAnalysis::InitializeCandidate( - const AnalysisRequest& request) { - // Clear dependent objects in reverse ownership order so a reused analysis - // never exposes a view into a Domain from an earlier run. - effective_rhs_.reset(); - partitioned_stiffness_.reset(); - full_stiffness_.reset(); - state_.reset(); - dofs_.reset(); - element_view_.clear(); - elements_.clear(); - model_.reset(); - domain_.reset(); - diagnostics_.clear(); - request_ = request; - - const auto parsed = AbaqusInputReader{}.Read(request_.input_path); - if (!parsed.HasValue()) { - return parsed.GetStatus(); - } - auto domain = AbaqusDomainMapper{}.Map(parsed.Value()); - if (!domain.HasValue()) { - return domain.GetStatus(); +Status LinearStaticAnalysis::Run(const AnalysisRequest& request) { + Result prepared = + PrepareAnalysisInput(request.input_path); + if (!prepared.HasValue()) { + return prepared.GetStatus(); } - domain_ = std::make_unique(std::move(domain.Value())); - diagnostics_ = domain_->Warnings(); - SortDiagnostics(diagnostics_); - return Status::Ok(); -} - -Status LinearStaticAnalysis::BuildAnalysisModel() { - auto model = AnalysisModel::Create(*domain_); - if (!model.HasValue()) { - return model.GetStatus(); - } - model_ = std::make_unique(std::move(model.Value())); - return Status::Ok(); -} - -Status LinearStaticAnalysis::BuildDofMapAndSparsePattern() { - std::vector> elements; - ElementView element_view; - elements.reserve(model_->ActiveElements().size()); - element_view.reserve(model_->ActiveElements().size()); - const ElementFactory factory; - for (const EntityIndex element_index : model_->ActiveElements()) { - auto element = factory.Create(domain_->Elements()[element_index], *domain_); - if (!element.HasValue()) { - return element.GetStatus(); - } - element_view.emplace_back(*element.Value()); - elements.push_back(std::move(element.Value())); + LinearStaticProcedureExecutor executor(parallel_for_, linear_solver_); + Result execution = + executor.Execute(prepared.Value().domain); + if (!execution.HasValue()) { + return execution.GetStatus(); } - DofManager dofs; - const Status dof_status = dofs.Build(*model_, element_view); - if (!dof_status.IsOk()) { - return dof_status; - } - - elements_ = std::move(elements); - element_view_ = std::move(element_view); - dofs_ = std::make_unique(std::move(dofs)); - state_ = std::make_unique( - AnalysisState::Create(*dofs_, {"Step-1", 0U})); - return Status::Ok(); -} - -Status LinearStaticAnalysis::AssembleAndPartitionStiffness() { - auto stiffness = - SparseAssembler::Assemble(element_view_, *dofs_, parallel_for_); - if (!stiffness.HasValue()) { - return stiffness.GetStatus(); - } - full_stiffness_ = - std::make_unique(std::move(stiffness.Value())); - - auto partitioned = - EssentialConstraintPolicy{}.Partition(*full_stiffness_, *dofs_); - if (!partitioned.HasValue()) { - return partitioned.GetStatus(); - } - partitioned_stiffness_ = - std::make_unique(std::move(partitioned.Value())); - return Status::Ok(); -} - -Status LinearStaticAnalysis::FactorizeFreeSystem() { - // This call intentionally precedes all load assembly in Run. - return linear_solver_.Factorize(partitioned_stiffness_->k_ff); -} - -Status LinearStaticAnalysis::AssembleLoadsAndEffectiveRhs() { - auto full_load = LoadAssembler::AssembleFullNodalLoad(*model_, *dofs_); - if (!full_load.HasValue()) { - return full_load.GetStatus(); - } - state_->ExternalForce() = std::move(full_load.Value()); - - auto rhs = LoadAssembler::EffectiveFreeRhs(state_->ExternalForce(), - partitioned_stiffness_->k_fc, - dofs_->PrescribedValues(), *dofs_); - if (!rhs.HasValue()) { - return rhs.GetStatus(); - } - effective_rhs_ = std::make_unique(std::move(rhs.Value())); - return Status::Ok(); -} - -Status LinearStaticAnalysis::SubstituteAndReconstruct() { - Vector free_displacement{dofs_->FreeDofCount()}; - const Status solve_status = - linear_solver_.Solve(*effective_rhs_, free_displacement); - if (!solve_status.IsOk()) { - return solve_status; - } - - state_->Displacement() = EssentialConstraintPolicy{}.ReconstructFull( - free_displacement, dofs_->PrescribedValues(), *dofs_); - return Status::Ok(); -} - -Status LinearStaticAnalysis::RecoverAndWrite() { - const Status recovery_status = ResultRecovery::Recover( - *model_, element_view_, *dofs_, *full_stiffness_, state_->Displacement(), - state_->ExternalForce(), *state_); - if (!recovery_status.IsOk()) { - return recovery_status; - } - return results_writer_.Write(request_.output_path, *domain_, *state_, - diagnostics_); + return results_writer_.Write(request.output_path, prepared.Value().domain, + execution.Value().state, + execution.Value().diagnostics); } } // namespace fesa diff --git a/src/fesa/analysis/linear_static_procedure.cpp b/src/fesa/analysis/linear_static_procedure.cpp new file mode 100644 index 0000000..cbe74df --- /dev/null +++ b/src/fesa/analysis/linear_static_procedure.cpp @@ -0,0 +1,105 @@ +#include "fesa/analysis/linear_static_procedure.h" + +#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.h" +#include "fesa/elements/element_factory.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/results/result_recovery.h" +#include "fesa/solvers/linear/linear_solver.h" + +namespace fesa { + +LinearStaticProcedureExecutor::LinearStaticProcedureExecutor( + const ParallelFor& parallel_for, LinearSolver& linear_solver) + : parallel_for_(parallel_for), linear_solver_(linear_solver) {} + +Result LinearStaticProcedureExecutor::Execute( + const Domain& domain) const { + Result model = AnalysisModel::Create(domain); + if (!model.HasValue()) { + return Result::Failure(model.GetStatus()); + } + + std::vector> elements; + ElementView element_view; + elements.reserve(model.Value().ActiveElements().size()); + element_view.reserve(model.Value().ActiveElements().size()); + const ElementFactory factory; + for (const EntityIndex element_index : model.Value().ActiveElements()) { + Result> element = + factory.Create(domain.Elements()[element_index], domain); + if (!element.HasValue()) { + return Result::Failure(element.GetStatus()); + } + element_view.emplace_back(*element.Value()); + elements.push_back(std::move(element.Value())); + } + + DofManager dofs; + const Status dof_status = dofs.Build(model.Value(), element_view); + if (!dof_status.IsOk()) { + return Result::Failure(dof_status); + } + + AnalysisState state = AnalysisState::Create(dofs, {"Step-1", 0U}); + Result stiffness = + SparseAssembler::Assemble(element_view, dofs, parallel_for_); + if (!stiffness.HasValue()) { + return Result::Failure(stiffness.GetStatus()); + } + + Result partitions = + EssentialConstraintPolicy{}.Partition(stiffness.Value(), dofs); + if (!partitions.HasValue()) { + return Result::Failure(partitions.GetStatus()); + } + + const Status factorized = linear_solver_.Factorize(partitions.Value().k_ff); + if (!factorized.IsOk()) { + return Result::Failure(factorized); + } + + Result external = + LoadAssembler::AssembleFullNodalLoad(model.Value(), dofs); + if (!external.HasValue()) { + return Result::Failure(external.GetStatus()); + } + state.ExternalForce() = std::move(external.Value()); + + Result rhs = LoadAssembler::EffectiveFreeRhs( + state.ExternalForce(), partitions.Value().k_fc, dofs.PrescribedValues(), + dofs); + if (!rhs.HasValue()) { + return Result::Failure(rhs.GetStatus()); + } + + Vector free_displacement{dofs.FreeDofCount()}; + const Status solved = linear_solver_.Solve(rhs.Value(), free_displacement); + if (!solved.IsOk()) { + return Result::Failure(solved); + } + state.Displacement() = EssentialConstraintPolicy{}.ReconstructFull( + free_displacement, dofs.PrescribedValues(), dofs); + + const Status recovered = ResultRecovery::Recover( + model.Value(), element_view, dofs, stiffness.Value(), + state.Displacement(), state.ExternalForce(), state); + if (!recovered.IsOk()) { + return Result::Failure(recovered); + } + + std::vector diagnostics = domain.Warnings(); + SortDiagnostics(diagnostics); + return Result::Success( + {std::move(state), std::move(diagnostics)}); +} + +} // namespace fesa diff --git a/src/fesa/analysis/nonlinear_static_procedure.cpp b/src/fesa/analysis/nonlinear_static_procedure.cpp new file mode 100644 index 0000000..41d79b8 --- /dev/null +++ b/src/fesa/analysis/nonlinear_static_procedure.cpp @@ -0,0 +1,807 @@ +#include "fesa/analysis/nonlinear_static_procedure.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/analysis/newton_iteration_log.h" +#include "fesa/analysis/nonlinear_static_control.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/constraints/essential_constraint_policy.h" +#include "fesa/elements/finite_rotation.h" +#include "fesa/elements/mitc4_shell.h" +#include "fesa/elements/nonlinear_mitc4_shell.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/math/sparse_matrix.h" +#include "fesa/properties/shell_section.h" +#include "fesa/solvers/linear/tangent_solver.h" + +namespace fesa { +namespace { + +constexpr std::size_t kDofsPerNode = 6U; +constexpr std::size_t kTranslationsPerNode = 3U; +constexpr std::size_t kShellDofCount = 24U; +constexpr std::uint32_t kMaximumNewtonIterations = 100U; + +struct RuntimeShell { + EntityIndex element_index; + std::array scatter; + NonlinearMitc4Shell element; +}; + +struct InternalAssembly { + Vector force; + SparseMatrix tangent; + std::vector shell_rows; + double physical_strain_energy; +}; + +Status ProcedureFailure(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}}); +} + +std::string DiagnosticCode(const Status& status) { + if (status.Diagnostics().empty()) { + return "nonlinear-procedure-failure"; + } + return status.Diagnostics().front().code; +} + +std::string DiagnosticReason(const Status& status) { + if (status.Diagnostics().empty()) { + return "The nonlinear procedure returned no diagnostic."; + } + return status.Diagnostics().front().message; +} + +NewtonChannelDecision UnevaluatedDecision(const std::string& identity, + const NewtonChannelUnit unit) { + return {identity, + 0.0, + 0.0, + NewtonConvergenceBranch::kZeroScale, + kNewtonZeroScaleThreshold, + unit, + false}; +} + +NewtonIterationLogRecord NewLogRecord(const StepDefinition& step, + const LoadFactorIncrement& increment, + const std::uint32_t iteration) { + return {step.Name(), + increment.increment_index, + increment.load_factor, + iteration, + NewtonIterationStage::kEntered, + {false, NewtonIterationStage::kEntered, {}, {}}, + "entered", + "pending", + {}, + UnevaluatedDecision("residual_force", NewtonChannelUnit::kForce), + UnevaluatedDecision("residual_moment", + NewtonChannelUnit::kForceLength), + UnevaluatedDecision("correction_translation", + NewtonChannelUnit::kLength), + UnevaluatedDecision("correction_rotation", + NewtonChannelUnit::kDimensionless), + false}; +} + +Status WriteAndFlush(NewtonIterationLogWriter& log, + const NewtonIterationLogRecord& record) { + const Status write_status = log.Write(record); + if (!write_status.IsOk()) { + return write_status; + } + return log.Flush(); +} + +Status RecordFailure(NewtonIterationLogWriter& log, + NewtonIterationLogRecord* record, const Status& failure, + const NewtonIterationStage stage) { + record->last_completed_stage = stage; + record->failure = {true, stage, DiagnosticCode(failure), + DiagnosticReason(failure)}; + record->transition = "failed"; + record->outcome = "failed"; + record->diagnostic_code = DiagnosticCode(failure); + const Status log_status = WriteAndFlush(log, *record); + return log_status.IsOk() ? failure : log_status; +} + +Matrix InitialFrame(const ShellNodeInitialFrame& frame) { + Matrix result{3U, 3U}; + for (std::size_t row = 0U; row < kTranslationsPerNode; ++row) { + result(row, 0U) = frame.tangent_a[row]; + result(row, 1U) = frame.tangent_b[row]; + result(row, 2U) = frame.director[row]; + } + return result; +} + +Result BuildDofManager(const Domain& domain) { + Result model = AnalysisModel::Create(domain); + if (!model.HasValue()) { + return Result::Failure(model.GetStatus()); + } + return DofManager::Create(model.Value()); +} + +Result BuildBaseline(const Domain& domain) { + if (domain.ShellNodeInitialFrames().size() != domain.Nodes().size()) { + return Result::Failure(ProcedureFailure( + FailureCategory::kModel, "invalid-nonlinear-state", + {domain.SourcePath(), 0U}, + "NODE", domain.SourceContentIdentity(), + "Nonlinear analysis requires one initial shell frame per node.")); + } + + NonlinearStateBaselineInput baseline; + baseline.node_baselines.reserve(domain.Nodes().size()); + for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { + const ShellNodeInitialFrame& frame = domain.ShellNodeInitialFrames()[node]; + if (frame.node_index != static_cast(node)) { + return Result::Failure(ProcedureFailure( + FailureCategory::kModel, "invalid-nonlinear-state", + domain.Nodes()[node].location, "NODE", + domain.Nodes()[node].source_id.source_label_text, + "Initial shell frames must follow stable Domain node order.")); + } + baseline.node_baselines.push_back( + {static_cast(node), domain.Nodes()[node].source_id, + {frame.tangent_a, frame.tangent_b, frame.director}, frame.director}); + } + return Result::Success(std::move(baseline)); +} + +Result> BuildRuntimeShells(const Domain& domain, + const DofManager& dofs) { + std::vector shells; + shells.reserve(domain.Elements().Size()); + for (std::size_t index = 0U; index < domain.Elements().Size(); ++index) { + const auto* definition = + dynamic_cast(&domain.Elements()[index]); + if (definition == nullptr) { + return Result>::Failure(ProcedureFailure( + FailureCategory::kModel, "unsupported-nonlinear-element", + {}, "ELEMENT", + domain.Elements()[index].SourceId().source_label_text, + "The nonlinear static procedure supports MITC4 shell elements " + "only.")); + } + if (definition->section_index >= domain.Properties().Size() || + definition->material_index >= domain.Materials().Size()) { + return Result>::Failure(ProcedureFailure( + FailureCategory::kModel, "invalid-nonlinear-static-data", + definition->location, "ELEMENT", + definition->source_id.source_label_text, + "The nonlinear shell element property or material index is " + "invalid.")); + } + + const auto* section = + dynamic_cast( + &domain.Properties()[definition->section_index]); + const auto* material = dynamic_cast( + &domain.Materials()[definition->material_index]); + if (section == nullptr || material == nullptr) { + return Result>::Failure(ProcedureFailure( + FailureCategory::kModel, "incompatible-nonlinear-shell-property", + definition->location, "ELEMENT", + definition->source_id.source_label_text, + "A nonlinear MITC4 shell requires a shell section and linear " + "elastic material.")); + } + + std::array nodes{}; + std::array frames{}; + for (std::size_t node = 0U; node < definition->node_indices.size(); + ++node) { + const EntityIndex node_index = definition->node_indices[node]; + if (node_index >= domain.Nodes().size() || + node_index >= domain.ShellNodeInitialFrames().size()) { + return Result>::Failure(ProcedureFailure( + FailureCategory::kModel, "invalid-shell-director", + definition->location, "ELEMENT", + definition->source_id.source_label_text, + "A nonlinear shell node does not have a stable initial frame.")); + } + nodes[node] = &domain.Nodes()[node_index]; + frames[node] = domain.ShellNodeInitialFrames()[node_index]; + if (frames[node].node_index != node_index) { + return Result>::Failure(ProcedureFailure( + FailureCategory::kModel, "invalid-shell-director", + definition->location, "ELEMENT", + definition->source_id.source_label_text, + "A nonlinear shell initial frame does not match its node " + "identity.")); + } + } + + Result shell = + NonlinearMitc4Shell::Create(nodes, frames, *section, *material); + if (!shell.HasValue()) { + return Result>::Failure(shell.GetStatus()); + } + const EntityIndex element_index = static_cast(index); + shells.push_back( + {element_index, dofs.ShellElementScatter(element_index), + std::move(shell.Value())}); + } + return Result>::Success(std::move(shells)); +} + +Status PopulateNodalKinematics(const Domain& domain, + NonlinearTrialState* trial) { + if (trial->generalized_coordinates.Size() != + domain.Nodes().size() * kDofsPerNode) { + return ProcedureFailure( + FailureCategory::kSolver, "nonfinite-newton-data", + {domain.SourcePath(), 0U}, + "NONLINEAR_STATIC_PROCEDURE", domain.SourceContentIdentity(), + "Nonlinear trial coordinates do not match the Domain node inventory."); + } + + trial->displacement = trial->generalized_coordinates; + trial->rotation_vector = Vector{trial->generalized_coordinates.Size()}; + trial->nodal_triad.clear(); + trial->nodal_director.clear(); + trial->nodal_triad.reserve(domain.Nodes().size()); + trial->nodal_director.reserve(domain.Nodes().size()); + for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { + const std::size_t first = node * kDofsPerNode; + const Vector3 phi{trial->generalized_coordinates[first + 3U], + trial->generalized_coordinates[first + 4U], + trial->generalized_coordinates[first + 5U]}; + Result rotation = MakeFiniteRotationState( + phi, InitialFrame(domain.ShellNodeInitialFrames()[node])); + if (!rotation.HasValue()) { + return rotation.GetStatus(); + } + for (std::size_t component = 0U; component < kTranslationsPerNode; + ++component) { + trial->rotation_vector[first + 3U + component] = phi[component]; + } + const Matrix& q = rotation.Value().q; + trial->nodal_triad.push_back( + {{{q(0U, 0U), q(1U, 0U), q(2U, 0U)}, + {q(0U, 1U), q(1U, 1U), q(2U, 1U)}, + {q(0U, 2U), q(1U, 2U), q(2U, 2U)}}}); + trial->nodal_director.push_back(rotation.Value().director.Components()); + } + return Status::Ok(); +} + +NonlinearMitc4TrialDofs ExtractElementTrial( + const RuntimeShell& shell, const Vector& coordinates) { + NonlinearMitc4TrialDofs trial{}; + for (std::size_t node = 0U; node < 4U; ++node) { + const std::size_t first = node * kDofsPerNode; + trial.displacement[node] = + {coordinates[shell.scatter[first]], + coordinates[shell.scatter[first + 1U]], + coordinates[shell.scatter[first + 2U]]}; + trial.rotation_vector[node] = + {coordinates[shell.scatter[first + 3U]], + coordinates[shell.scatter[first + 4U]], + coordinates[shell.scatter[first + 5U]]}; + } + return trial; +} + +bool IsFinite(const Vector& values) { + for (std::size_t index = 0U; index < values.Size(); ++index) { + if (!std::isfinite(values[index])) { + return false; + } + } + return true; +} + +bool IsFinite(const Matrix& values) { + for (std::size_t row = 0U; row < values.Rows(); ++row) { + for (std::size_t column = 0U; column < values.Columns(); ++column) { + if (!std::isfinite(values(row, column))) { + return false; + } + } + } + return true; +} + +Result AssembleInternal( + const std::vector& shells, const DofManager& dofs, + const Vector& coordinates, const ParallelFor& parallel_for) { + std::vector> local( + shells.size()); + std::vector> failures(shells.size()); + parallel_for.Execute(shells.size(), [&](const std::size_t index) { + Result result = + EvaluateInternalResidualAndTangent(shells[index].element, + ExtractElementTrial(shells[index], + coordinates)); + if (!result.HasValue()) { + failures[index] = result.GetStatus(); + return; + } + local[index] = std::move(result.Value()); + }); + + Vector force{dofs.FullDofCount()}; + std::vector contributions; + contributions.reserve(shells.size() * kShellDofCount * kShellDofCount); + std::vector rows; + rows.reserve(shells.size() * 4U); + double energy = 0.0; + for (std::size_t index = 0U; index < shells.size(); ++index) { + if (failures[index].has_value()) { + return Result::Failure(*failures[index]); + } + if (!local[index].has_value() || + !IsFinite(local[index]->internal_residual) || + !IsFinite(local[index]->internal_tangent) || + !std::isfinite(local[index]->physical_strain_energy)) { + return Result::Failure(ProcedureFailure( + FailureCategory::kSolver, "nonfinite-newton-data", {}, + "NONLINEAR_STATIC_PROCEDURE", std::to_string(index), + "A nonlinear shell internal residual or tangent is nonfinite.")); + } + for (std::size_t row = 0U; row < kShellDofCount; ++row) { + force[shells[index].scatter[row]] += local[index]->internal_residual[row]; + for (std::size_t column = 0U; column < kShellDofCount; ++column) { + contributions.push_back( + {shells[index].scatter[row], shells[index].scatter[column], + local[index]->internal_tangent(row, column), index, + row * kShellDofCount + column}); + } + } + for (const NonlinearMitc4RecoveryPoint& point : + local[index]->recovery.gauss_points) { + rows.push_back({shells[index].element_index, + point.location, + point.natural_coordinates, + {point.local_frame.e1, point.local_frame.e2, + point.local_frame.e3}, + point.generalized_strain, + point.section_resultant, + point.stress}); + } + energy += local[index]->physical_strain_energy; + } + if (!IsFinite(force) || !std::isfinite(energy)) { + return Result::Failure(ProcedureFailure( + FailureCategory::kSolver, "nonfinite-newton-data", {}, + "NONLINEAR_STATIC_PROCEDURE", {}, + "Global nonlinear internal force or energy is nonfinite.")); + } + + Result 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(rows), energy}); +} + +Result Difference(const Vector& left, const Vector& right) { + if (left.Size() != right.Size()) { + return Result::Failure(ProcedureFailure( + FailureCategory::kSolver, "nonfinite-newton-data", {}, + "NONLINEAR_STATIC_PROCEDURE", {}, + "Nonlinear internal and external vectors must have equal dimensions.")); + } + Vector result{left.Size()}; + for (std::size_t index = 0U; index < result.Size(); ++index) { + result[index] = left[index] - right[index]; + } + if (!IsFinite(result)) { + return Result::Failure(ProcedureFailure( + FailureCategory::kSolver, "nonfinite-newton-data", {}, + "NONLINEAR_STATIC_PROCEDURE", {}, + "Nonlinear residual must remain finite.")); + } + return Result::Success(std::move(result)); +} + +Result Difference(const SparseMatrix& internal, + const SparseMatrix& external, + const DofManager& dofs) { + if (internal.Rows() != external.Rows() || + internal.Columns() != external.Columns() || + internal.RowOffsets() != external.RowOffsets() || + internal.ColumnIndices() != external.ColumnIndices()) { + return Result::Failure(ProcedureFailure( + FailureCategory::kSolver, "nonfinite-newton-data", {}, + "NONLINEAR_STATIC_PROCEDURE", {}, + "Nonlinear internal and external tangent patterns must match.")); + } + std::vector contributions; + contributions.reserve(internal.Values().size()); + for (std::size_t row = 0U; row < internal.Rows(); ++row) { + for (std::size_t offset = internal.RowOffsets()[row]; + offset < internal.RowOffsets()[row + 1U]; ++offset) { + const double value = + internal.Values()[offset] - external.Values()[offset]; + if (!std::isfinite(value)) { + return Result::Failure(ProcedureFailure( + FailureCategory::kSolver, "nonfinite-newton-data", {}, + "NONLINEAR_STATIC_PROCEDURE", {}, + "Nonlinear tangent must remain finite.")); + } + contributions.push_back( + {row, internal.ColumnIndices()[offset], value, row, offset}); + } + } + return SparseMatrix::FromCoo(internal.Rows(), internal.Columns(), + std::move(contributions), + dofs.GetSparsePattern()); +} + +Vector SelectFreeComponents(const Vector& full, const DofManager& dofs, + const bool rotations) { + std::size_t count = 0U; + for (const std::size_t full_dof : dofs.FreeDofs()) { + if ((full_dof % kDofsPerNode >= kTranslationsPerNode) == rotations) { + ++count; + } + } + Vector selected{count}; + std::size_t output = 0U; + for (const std::size_t full_dof : dofs.FreeDofs()) { + if ((full_dof % kDofsPerNode >= kTranslationsPerNode) == rotations) { + selected[output++] = full[full_dof]; + } + } + return selected; +} + +Vector SelectCorrectionComponents(const Vector& correction, + const DofManager& dofs, + const bool rotations) { + std::size_t count = 0U; + for (const std::size_t full_dof : dofs.FreeDofs()) { + if ((full_dof % kDofsPerNode >= kTranslationsPerNode) == rotations) { + ++count; + } + } + Vector selected{count}; + std::size_t output = 0U; + for (std::size_t equation = 0U; equation < dofs.FreeDofs().size(); + ++equation) { + if ((dofs.FreeDofs()[equation] % kDofsPerNode >= + kTranslationsPerNode) == rotations) { + selected[output++] = correction[equation]; + } + } + return selected; +} + +Vector CurrentIncrementComponents(const Vector& current, const Vector& start, + const DofManager& dofs, + const bool rotations) { + Vector difference{current.Size()}; + for (std::size_t index = 0U; index < difference.Size(); ++index) { + difference[index] = current[index] - start[index]; + } + return SelectFreeComponents(difference, dofs, rotations); +} + +std::array GlobalEquilibrium(const Domain& domain, + const Vector& residual) { + std::array equilibrium{}; + for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { + const std::size_t first = node * kDofsPerNode; + const Vector3 force{residual[first], residual[first + 1U], + residual[first + 2U]}; + const Vector3 position{domain.Nodes()[node].coordinates}; + const Vector3 moment = position.Cross(force); + equilibrium[0U] += force.X(); + equilibrium[1U] += force.Y(); + equilibrium[2U] += force.Z(); + equilibrium[3U] += moment.X() + residual[first + 3U]; + equilibrium[4U] += moment.Y() + residual[first + 4U]; + equilibrium[5U] += moment.Z() + residual[first + 5U]; + } + return equilibrium; +} + +Status ApplyIncrementConstraints(const StepDefinition& step, + const LoadFactorIncrement& increment, + const DofManager& dofs, + NonlinearTrialState* trial) { + Result source_targets = + ApplyAbsolutePrescribedDisplacementTarget(step, increment.load_factor); + if (!source_targets.HasValue()) { + return source_targets.GetStatus(); + } + for (std::size_t index = 0U; index < dofs.ConstrainedDofs().size(); + ++index) { + const double target = + increment.load_factor * dofs.PrescribedValues()[index]; + if (!std::isfinite(target)) { + return ProcedureFailure( + FailureCategory::kSolver, "nonfinite-newton-data", step.Location(), + "BOUNDARY", step.Name(), + "A nonlinear prescribed displacement target is nonfinite."); + } + trial->generalized_coordinates[dofs.ConstrainedDofs()[index]] = target; + } + return Status::Ok(); +} + +AcceptedFrameState AcceptedFrame(const NonlinearAnalysisState& state, + const NonlinearTrialState& trial, + const LoadFactorIncrement& increment, + const std::uint32_t iteration) { + return {static_cast(state.AcceptedFrames().size()), + increment.increment_index, + increment.load_factor, + iteration, + increment.is_factor_one, + trial.nodal_internal_order, + trial.nodal_source_ids, + trial.generalized_coordinates, + trial.displacement, + trial.rotation_vector, + trial.nodal_triad, + trial.nodal_director, + trial.internal_force, + trial.external_force, + trial.residual, + trial.reaction, + trial.shell_results, + trial.physical_strain_energy, + trial.global_equilibrium}; +} + +} // namespace + +NonlinearStaticProcedureExecutor::NonlinearStaticProcedureExecutor( + const ParallelFor& parallel_for, TangentSolver& tangent_solver) + : parallel_for_(parallel_for), tangent_solver_(tangent_solver) {} + +Result +NonlinearStaticProcedureExecutor::Execute( + const Domain& domain, NewtonIterationLogWriter& log) const { + if (domain.Steps().Size() != 1U || !domain.Steps()[0U].IsNonlinearStatic()) { + return Result::Failure(ProcedureFailure( + FailureCategory::kInput, "invalid-nonlinear-static-data", + {domain.SourcePath(), 0U}, "STATIC", domain.SourceContentIdentity(), + "Nonlinear procedure execution requires exactly one NLGEOM static " + "step.")); + } + const StepDefinition& step = domain.Steps()[0U]; + + Result dofs = BuildDofManager(domain); + if (!dofs.HasValue()) { + return Result::Failure(dofs.GetStatus()); + } + Result baseline = BuildBaseline(domain); + if (!baseline.HasValue()) { + return Result::Failure( + baseline.GetStatus()); + } + Result state = + NonlinearAnalysisState::Create(domain, dofs.Value(), baseline.Value()); + if (!state.HasValue()) { + return Result::Failure(state.GetStatus()); + } + Result> shells = + BuildRuntimeShells(domain, dofs.Value()); + if (!shells.HasValue()) { + return Result::Failure(shells.GetStatus()); + } + Result> increments = + BuildFixedLoadFactorSequence(step); + if (!increments.HasValue()) { + return Result::Failure( + increments.GetStatus()); + } + + const EssentialConstraintPolicy constraints; + for (const LoadFactorIncrement& increment : increments.Value()) { + Result trial = state.Value().BeginTrialFromCommitted(); + if (!trial.HasValue()) { + return Result::Failure(trial.GetStatus()); + } + const Vector increment_start = trial.Value().generalized_coordinates; + Status constraint_status = + ApplyIncrementConstraints(step, increment, dofs.Value(), + &trial.Value()); + if (!constraint_status.IsOk()) { + return Result::Failure(constraint_status); + } + + for (std::uint32_t iteration = 1U; iteration <= kMaximumNewtonIterations; + ++iteration) { + NewtonIterationLogRecord record = + NewLogRecord(step, increment, iteration); + trial.Value().increment_index = increment.increment_index; + trial.Value().load_factor = increment.load_factor; + trial.Value().is_factor_zero_baseline = false; + + Status kinematics = PopulateNodalKinematics(domain, &trial.Value()); + if (!kinematics.IsOk()) { + return Result::Failure( + RecordFailure(log, &record, kinematics, + NewtonIterationStage::kKinematicsEvaluated)); + } + record.last_completed_stage = NewtonIterationStage::kKinematicsEvaluated; + + Result internal = AssembleInternal( + shells.Value(), dofs.Value(), trial.Value().generalized_coordinates, + parallel_for_); + if (!internal.HasValue()) { + return Result::Failure(RecordFailure( + log, &record, internal.GetStatus(), + NewtonIterationStage::kInternalAssembled)); + } + record.last_completed_stage = NewtonIterationStage::kInternalAssembled; + + Result external = + AssembleGlobalFixedLoadContribution( + domain, dofs.Value(), trial.Value().generalized_coordinates, + increment.load_factor); + if (!external.HasValue()) { + return Result::Failure(RecordFailure( + log, &record, external.GetStatus(), + NewtonIterationStage::kExternalAssembled)); + } + record.last_completed_stage = NewtonIterationStage::kExternalAssembled; + + Result residual = + Difference(internal.Value().force, + external.Value().generalized_external_force); + if (!residual.HasValue()) { + return Result::Failure(RecordFailure( + log, &record, residual.GetStatus(), + NewtonIterationStage::kExternalAssembled)); + } + Result tangent = + Difference(internal.Value().tangent, + external.Value().generalized_external_tangent, + dofs.Value()); + if (!tangent.HasValue()) { + return Result::Failure(RecordFailure( + log, &record, tangent.GetStatus(), + NewtonIterationStage::kExternalAssembled)); + } + trial.Value().internal_force = std::move(internal.Value().force); + trial.Value().external_force = + std::move(external.Value().generalized_external_force); + trial.Value().residual = std::move(residual.Value()); + trial.Value().reaction = trial.Value().residual; + trial.Value().shell_results = std::move(internal.Value().shell_rows); + trial.Value().physical_strain_energy = + internal.Value().physical_strain_energy; + trial.Value().global_equilibrium = + GlobalEquilibrium(domain, trial.Value().residual); + + Result partition = + constraints.Partition(tangent.Value(), dofs.Value()); + if (!partition.HasValue()) { + return Result::Failure(RecordFailure( + log, &record, partition.GetStatus(), + NewtonIterationStage::kExternalAssembled)); + } + const Status factorization = tangent_solver_.FactorizeFresh( + partition.Value().k_ff); + if (!factorization.IsOk()) { + return Result::Failure(RecordFailure( + log, &record, factorization, + NewtonIterationStage::kSolvedCorrection)); + } + + Vector correction_rhs = constraints.GatherFree(trial.Value().residual, + dofs.Value()); + for (std::size_t index = 0U; index < correction_rhs.Size(); ++index) { + correction_rhs[index] = -correction_rhs[index]; + } + Vector correction{dofs.Value().FreeDofCount()}; + const Status solve = + tangent_solver_.SolveCorrection(correction_rhs, correction); + if (!solve.IsOk()) { + return Result::Failure(RecordFailure( + log, &record, solve, NewtonIterationStage::kSolvedCorrection)); + } + record.last_completed_stage = NewtonIterationStage::kSolvedCorrection; + + const NewtonConvergenceResult convergence = EvaluateNewtonConvergence({ + {"residual_force", SelectFreeComponents(trial.Value().residual, + dofs.Value(), false), + SelectFreeComponents(trial.Value().internal_force, dofs.Value(), + false), + SelectFreeComponents(trial.Value().external_force, dofs.Value(), + false)}, + {"residual_moment", SelectFreeComponents(trial.Value().residual, + dofs.Value(), true), + SelectFreeComponents(trial.Value().internal_force, dofs.Value(), + true), + SelectFreeComponents(trial.Value().external_force, dofs.Value(), + true)}, + {"correction_translation", + SelectCorrectionComponents(correction, dofs.Value(), false), + CurrentIncrementComponents(trial.Value().generalized_coordinates, + increment_start, dofs.Value(), false)}, + {"correction_rotation", + SelectCorrectionComponents(correction, dofs.Value(), true), + CurrentIncrementComponents(trial.Value().generalized_coordinates, + increment_start, dofs.Value(), true)}}); + record.force_residual = convergence.force_residual; + record.moment_residual = convergence.moment_residual; + record.translation_correction = convergence.translation_correction; + record.rotation_correction = convergence.rotation_correction; + record.combined_passed = convergence.combined_passed; + record.last_completed_stage = NewtonIterationStage::kConvergenceEvaluated; + + if (convergence.combined_passed) { + const Status committed = state.Value().CommitAcceptedFrame( + AcceptedFrame(state.Value(), trial.Value(), increment, iteration)); + if (!committed.IsOk()) { + return Result::Failure(RecordFailure( + log, &record, committed, NewtonIterationStage::kAccepted)); + } + record.last_completed_stage = NewtonIterationStage::kAccepted; + record.transition = "accepted"; + record.outcome = "accepted"; + const Status log_status = WriteAndFlush(log, record); + if (!log_status.IsOk()) { + return Result::Failure(log_status); + } + break; + } + + if (iteration == kMaximumNewtonIterations) { + const Status failure = ProcedureFailure( + FailureCategory::kSolver, "nonlinear-nonconvergence", + step.Location(), "STATIC", step.Name(), + "The nonlinear procedure did not converge within 100 iterations."); + return Result::Failure(RecordFailure( + log, &record, failure, + NewtonIterationStage::kConvergenceEvaluated)); + } + + for (std::size_t equation = 0U; equation < dofs.Value().FreeDofs().size(); + ++equation) { + const std::size_t full_dof = dofs.Value().FreeDofs()[equation]; + trial.Value().generalized_coordinates[full_dof] += correction[equation]; + if (!std::isfinite(trial.Value().generalized_coordinates[full_dof])) { + const Status failure = ProcedureFailure( + FailureCategory::kSolver, "nonfinite-newton-data", {}, + "NONLINEAR_STATIC_PROCEDURE", std::to_string(full_dof), + "A nonlinear free-coordinate correction is nonfinite."); + return Result::Failure(RecordFailure( + log, &record, failure, NewtonIterationStage::kSolvedCorrection)); + } + } + record.transition = "continue"; + record.outcome = "not-converged"; + const Status log_status = WriteAndFlush(log, record); + if (!log_status.IsOk()) { + return Result::Failure(log_status); + } + } + } + + std::vector diagnostics = domain.Warnings(); + SortDiagnostics(diagnostics); + return Result::Success( + {std::move(state.Value()), std::move(diagnostics)}); +} + +} // namespace fesa diff --git a/src/fesa/analysis/prepared_analysis.cpp b/src/fesa/analysis/prepared_analysis.cpp new file mode 100644 index 0000000..ce2042d --- /dev/null +++ b/src/fesa/analysis/prepared_analysis.cpp @@ -0,0 +1,79 @@ +#include "fesa/analysis/prepared_analysis.h" + +#include + +#include "fesa/analysis/linear_static_procedure.h" +#include "fesa/analysis/newton_iteration_log.h" +#include "fesa/analysis/nonlinear_static_procedure.h" +#include "fesa/app/output_paths.h" +#include "fesa/io/abaqus/domain_mapper.h" +#include "fesa/io/abaqus/input_reader.h" +#include "fesa/results/nonlinear_results_writer.h" +#include "fesa/results/results_writer.h" + +namespace fesa { + +Result PrepareAnalysisInput( + const std::filesystem::path& input_path) { + const Result syntax = AbaqusInputReader{}.Read(input_path); + if (!syntax.HasValue()) { + return Result::Failure(syntax.GetStatus()); + } + + Result mapped = AbaqusDomainMapper{}.Map(syntax.Value()); + if (!mapped.HasValue()) { + return Result::Failure(mapped.GetStatus()); + } + + PreparedAnalysisInput prepared{std::move(mapped.Value()), {}}; + prepared.diagnostics = prepared.domain.Warnings(); + SortDiagnostics(prepared.diagnostics); + return Result::Success(std::move(prepared)); +} + +Result ExecutePreparedAnalysis( + const Domain& domain, LinearStaticProcedureExecutor& linear_executor, + NonlinearStaticProcedureExecutor& nonlinear_executor, + NewtonIterationLogWriter& log) { + if (domain.Steps().Size() != 1U) { + return Result::Failure(Status::Failure( + FailureCategory::kInput, + {{Severity::kError, "invalid-analysis-step-count", {}, "*STEP", {}, + "Prepared analysis requires exactly one static step."}})); + } + + if (domain.Steps()[0U].IsNonlinearStatic()) { + Result nonlinear = + nonlinear_executor.Execute(domain, log); + if (!nonlinear.HasValue()) { + return Result::Failure(nonlinear.GetStatus()); + } + return Result::Success( + PreparedAnalysisOutcome{std::move(nonlinear.Value())}); + } + + Result linear = linear_executor.Execute(domain); + if (!linear.HasValue()) { + return Result::Failure(linear.GetStatus()); + } + return Result::Success( + PreparedAnalysisOutcome{std::move(linear.Value())}); +} + +Status FinalizePreparedAnalysisResult( + const DerivedOutputPaths& paths, const Domain& domain, + const PreparedAnalysisOutcome& outcome, + const std::vector& diagnostics, ResultsWriter& linear_writer, + NonlinearResultsWriter& nonlinear_writer) { + if (std::holds_alternative(outcome)) { + return linear_writer.Write( + paths.results_path, domain, + std::get(outcome).state, diagnostics); + } + + return nonlinear_writer.Write( + paths.results_path, domain, + std::get(outcome).state, diagnostics); +} + +} // namespace fesa diff --git a/src/fesa/app/fesa_application.cpp b/src/fesa/app/fesa_application.cpp index e878035..5ac723b 100644 --- a/src/fesa/app/fesa_application.cpp +++ b/src/fesa/app/fesa_application.cpp @@ -2,17 +2,20 @@ #include #include -#include #include #include -#include "fesa/analysis/analysis.h" -#include "fesa/analysis/linear_static_analysis.h" +#include "fesa/analysis/linear_static_procedure.h" +#include "fesa/analysis/newton_iteration_log.h" +#include "fesa/analysis/nonlinear_static_procedure.h" +#include "fesa/analysis/prepared_analysis.h" #include "fesa/app/output_paths.h" #include "fesa/assembly/parallel_for.h" #include "fesa/core/diagnostic.h" #include "fesa/io/hdf5/hdf5_results_writer.h" +#include "fesa/io/hdf5/nonlinear_hdf5_results_writer.h" #include "fesa/solvers/linear/mkl_pardiso_solver.h" +#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" namespace fesa { namespace { @@ -73,17 +76,53 @@ int FesaApplication::Run(const std::vector& arguments) { return kOutputExitCode; } - AnalysisRequest request; - request.input_path = paths.Value().input_path; - request.output_path = paths.Value().results_path; - request.log_path = paths.Value().log_path; - TbbParallelFor parallel_for; MklPardisoSolver linear_solver; - Hdf5ResultsWriter results_writer; - std::unique_ptr analysis = std::make_unique( - parallel_for, linear_solver, results_writer); - const Status status = analysis->Run(request); + MklPardisoTangentSolver tangent_solver; + Hdf5ResultsWriter linear_writer; + NonlinearHdf5ResultsWriter nonlinear_writer; + NewtonIterationLogWriter log; + + Status status = log.Open(paths.Value().log_path); + if (!status.IsOk()) { + WriteDiagnostics(status.Diagnostics()); + return ExitCodeFor(status); + } + + Result prepared = + PrepareAnalysisInput(paths.Value().input_path); + if (!prepared.HasValue()) { + status = prepared.GetStatus(); + } else { + LinearStaticProcedureExecutor linear_executor(parallel_for, linear_solver); + NonlinearStaticProcedureExecutor nonlinear_executor(parallel_for, + tangent_solver); + Result outcome = ExecutePreparedAnalysis( + prepared.Value().domain, linear_executor, nonlinear_executor, log); + if (!outcome.HasValue()) { + status = outcome.GetStatus(); + } else { + const Status close_status = log.Close(); + if (!close_status.IsOk()) { + WriteDiagnostics(close_status.Diagnostics()); + return ExitCodeFor(close_status); + } + status = FinalizePreparedAnalysisResult( + paths.Value(), prepared.Value().domain, outcome.Value(), + prepared.Value().diagnostics, linear_writer, nonlinear_writer); + if (status.IsOk()) { + return kSuccessExitCode; + } + WriteDiagnostics(status.Diagnostics()); + return ExitCodeFor(status); + } + } + + const Status close_status = log.Close(); + if (!close_status.IsOk()) { + WriteDiagnostics(close_status.Diagnostics()); + return ExitCodeFor(close_status); + } if (status.IsOk()) { return kSuccessExitCode; } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 43fc174..3298c37 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -72,6 +72,7 @@ target_link_libraries( add_executable( fesa_integration_tests integration/analysis/linear_static_analysis_test.cpp + integration/analysis/nonlinear_static_procedure_test.cpp integration/app/fesa_application_test.cpp ) diff --git a/tests/integration/analysis/linear_static_analysis_test.cpp b/tests/integration/analysis/linear_static_analysis_test.cpp index a738210..2205230 100644 --- a/tests/integration/analysis/linear_static_analysis_test.cpp +++ b/tests/integration/analysis/linear_static_analysis_test.cpp @@ -16,6 +16,7 @@ #include "fesa/analysis/analysis.h" #include "fesa/assembly/parallel_for.h" +#include "fesa/math/sparse_matrix.h" #include "fesa/results/results_writer.h" #include "fesa/solvers/linear/mkl_pardiso_solver.h" diff --git a/tests/integration/analysis/nonlinear_static_procedure_test.cpp b/tests/integration/analysis/nonlinear_static_procedure_test.cpp new file mode 100644 index 0000000..dd13cab --- /dev/null +++ b/tests/integration/analysis/nonlinear_static_procedure_test.cpp @@ -0,0 +1,157 @@ +#include "fesa/app/fesa_application.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace { + +class TempDirectory { + public: + explicit TempDirectory(const std::string& label) { + static std::atomic sequence{0U}; + const auto tick = + std::chrono::steady_clock::now().time_since_epoch().count(); + path_ = std::filesystem::temp_directory_path() / + ("fesa-step10-nonlinear-" + label + "-" + + std::to_string(tick) + "-" + + std::to_string(sequence.fetch_add(1U))); + std::error_code error; + if (!std::filesystem::create_directory(path_, error) || error) { + throw std::runtime_error{"Unable to create nonlinear procedure fixture."}; + } + } + + TempDirectory(const TempDirectory&) = delete; + TempDirectory& operator=(const TempDirectory&) = delete; + + ~TempDirectory() { + std::error_code ignored; + std::filesystem::remove_all(path_, ignored); + } + + const std::filesystem::path& Path() const noexcept { return path_; } + + private: + std::filesystem::path path_; +}; + +void WriteText(const std::filesystem::path& path, const std::string& text) { + std::ofstream stream{path, std::ios::binary | std::ios::trunc}; + stream.write(text.data(), static_cast(text.size())); + if (!stream) { + throw std::runtime_error{"Unable to write nonlinear procedure input."}; + } +} + +std::string AllConstrainedNonlinearShellDeck() { + return R"inp(*Part, name=ShellPart +*Node +1, 0., 0., 0. +2, 1., 0., 0. +3, 1., 1., 0. +4, 0., 1., 0. +*Element, type=S4 +1, 1, 2, 3, 4 +*Elset, elset=ShellSet +1 +*Shell Section, elset=ShellSet, material=Steel +0.1 +*End Part +*Assembly, name=Assembly +*Instance, name=Shell-1, part=ShellPart +*End Instance +*Nset, nset=All, instance=Shell-1 +1, 2, 3, 4 +*End Assembly +*Material, name=Steel +*Elastic +1000., 0.25 +*Boundary +All, 1, 6 +*Step, name=Step-1, nlgeom=YES +*Static +0.1, 1., 0.001, 0.1 +*End Step +)inp"; +} + +class Hdf5File { + public: + explicit Hdf5File(const std::filesystem::path& path) + : handle_{H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)} { + if (handle_ < 0) { + throw std::runtime_error{"Unable to open nonlinear results HDF5."}; + } + } + + Hdf5File(const Hdf5File&) = delete; + Hdf5File& operator=(const Hdf5File&) = delete; + + ~Hdf5File() { (void)H5Fclose(handle_); } + + hid_t Get() const noexcept { return handle_; } + + private: + hid_t handle_; +}; + +std::string ReadStringAttribute(const hid_t object, const char* name) { + const hid_t attribute = H5Aopen(object, name, H5P_DEFAULT); + const hid_t type = attribute < 0 ? -1 : H5Aget_type(attribute); + if (attribute < 0 || type < 0) { + if (type >= 0) { + (void)H5Tclose(type); + } + if (attribute >= 0) { + (void)H5Aclose(attribute); + } + throw std::runtime_error{"Unable to read nonlinear metadata attribute."}; + } + char* value = nullptr; + const herr_t read = H5Aread(attribute, type, &value); + const std::string result = read >= 0 && value != nullptr ? value : ""; + if (value != nullptr) { + (void)H5free_memory(value); + } + (void)H5Tclose(type); + (void)H5Aclose(attribute); + if (read < 0) { + throw std::runtime_error{"Unable to read nonlinear metadata text."}; + } + return result; +} + +// NSNR-PROC-004, NSNR-PROC-007, and NSNR-PROC-010. +TEST(NonlinearStaticProcedure, + AllConstrainedNlgeomRunWritesSchemaV1AfterAcceptedFrames) { + TempDirectory directory{"schema-v1"}; + const auto input = directory.Path() / "shellnl.inp"; + const auto output = directory.Path() / "shellnl_results.h5"; + WriteText(input, AllConstrainedNonlinearShellDeck()); + + testing::internal::CaptureStderr(); + const int exit_code = fesa::FesaApplication{}.Run({input.string()}); + const std::string standard_error = testing::internal::GetCapturedStderr(); + + ASSERT_EQ(exit_code, 0) << standard_error; + ASSERT_TRUE(std::filesystem::exists(output)); + ASSERT_GT(H5Fis_hdf5(output.string().c_str()), 0); + const Hdf5File file{output}; + const hid_t metadata = H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT); + ASSERT_GE(metadata, 0); + EXPECT_EQ(ReadStringAttribute(metadata, "feature_id"), + "nonlinear-static-newton-raphson"); + (void)H5Gclose(metadata); + EXPECT_GT(H5Lexists(file.Get(), "/steps/Step-1/accepted_frames", H5P_DEFAULT), + 0); +} + +} // namespace diff --git a/tests/integration/app/fesa_application_test.cpp b/tests/integration/app/fesa_application_test.cpp index 330047c..efdbb4f 100644 --- a/tests/integration/app/fesa_application_test.cpp +++ b/tests/integration/app/fesa_application_test.cpp @@ -368,6 +368,20 @@ TEST(LinearStaticCli, DefaultOutputUsesInputSiblingDerivedName) { EXPECT_FALSE(std::filesystem::exists(directory.Path() / "results.h5")); } +// NSNR-PROC-001: log lifecycle begins after leasing and before parse/map. +TEST(FesaApplication, OpensDerivedLogBeforeInputMapping) { + TempDirectory directory{"log-before-input"}; + const auto input = directory.Path() / "missing.inp"; + const auto log = directory.Path() / "missing_log.txt"; + + const auto result = RunApplication({input.string()}); + + EXPECT_EQ(result.exit_code, 3); + ExpectDiagnosticFieldOrder(result.standard_error); + ASSERT_TRUE(std::filesystem::exists(log)); + EXPECT_EQ(std::filesystem::file_size(log), 0U); +} + TEST(LinearStaticCli, ReturnsEveryExactExitCodeAndOrderedDiagnostic) { TempDirectory directory{"exit-codes"}; const auto valid_input = directory.Path() / "valid.inp"; diff --git a/tests/reference/b33_reference_comparison_test.cpp b/tests/reference/b33_reference_comparison_test.cpp index 0c277f1..2987869 100644 --- a/tests/reference/b33_reference_comparison_test.cpp +++ b/tests/reference/b33_reference_comparison_test.cpp @@ -252,6 +252,7 @@ TEST(B33ReferenceComparison, GeneratesAuthoritativeHdf5AndComparisonEvidence) { std::sort(generated_names.begin(), generated_names.end()); EXPECT_EQ(generated_names, (std::vector{"cantilever beam.inp", + "cantilever beam_log.txt", "cantilever beam_results.h5", "comparison.json"})); diff --git a/tests/reference/mitc4_reference_cases_test.cpp b/tests/reference/mitc4_reference_cases_test.cpp index 360a070..75b3fe6 100644 --- a/tests/reference/mitc4_reference_cases_test.cpp +++ b/tests/reference/mitc4_reference_cases_test.cpp @@ -111,6 +111,7 @@ CaseEvidence RunCase(const std::string& case_id, EXPECT_EQ(generated, (std::vector{"comparison.json", input.filename().string(), + input.stem().string() + "_log.txt", results.filename().string()})); EXPECT_TRUE(std::filesystem::is_directory(reference_directory)); ExpectUnchanged(input, input_before);