FESADev/src/fesa/analysis/hht_dynamic_internal.h

#ifndef FESA_ANALYSIS_HHT_DYNAMIC_INTERNAL_H_
#define FESA_ANALYSIS_HHT_DYNAMIC_INTERNAL_H_

#include <array>
#include <cstdint>
#include <filesystem>
#include <memory>
#include <optional>
#include <vector>

#include "fesa/analysis/dynamic_analysis_state.h"
#include "fesa/analysis/hht_time_integration.h"
#include "fesa/elements/element_capabilities.h"
#include "fesa/elements/nonlinear_mitc4_inertia.h"
#include "fesa/elements/nonlinear_mitc4_shell.h"
#include "fesa/math/sparse_matrix.h"
#include "elements/nonlinear_mitc4_configuration_internal.h"
#include "math/dense_blas_internal.h"

namespace fesa {
class LinearSolver;
class ParallelFor;
class TangentSolver;
namespace hht_internal {

/// @brief Keeps immutable local evaluators with the DOF-owner scatter.
struct RuntimeShell {
  EntityIndex element_index;
  std::array<std::size_t, 24> scatter;
  std::unique_ptr<NonlinearMitc4Shell> internal_kernel;
  std::unique_ptr<NonlinearMitc4Inertia> inertia_kernel;
  std::unique_ptr<ShellTrialElement> internal;
  std::unique_ptr<ShellInertiaElement> inertia;
};

/// @brief Captures only controls used by the current MSVC x64 local evaluator.
struct TrialEvaluationContext {
  int rounding{-1};
  unsigned int mxcsr_control{0U};
  dense_blas_internal::LocalDensePolicy dense_policy{0, -1};
  bool supported{false};
};

/// @brief Queries metadata without changing FP controls or the restored TLS.
TrialEvaluationContext ReadTrialEvaluationContext() noexcept;

/// @brief Allows only masked, stable, single-worker nondynamic local policy.
bool ReusableTrialContext(const TrialEvaluationContext& before,
                          const TrialEvaluationContext& after) noexcept;

/// @brief Owns only complete directional constraint data for one local trial.
struct DirectionalShellConstraintData {
  EntityIndex element_index{0U};
  std::array<std::size_t, 24> scatter{};
  std::array<SourceEntityId, 4> node_source_ids{};
  nonlinear_mitc4_internal::ConfigurationSnapshot configuration;
  std::array<double, 24> q{};
  std::array<double, 24> v{};
  TrialEvaluationContext evaluation_context;
  bool directional_internal{false};
  Vector drilling_constraint{0U};
  Matrix drilling_constraint_gradient{0U, 0U};
  std::array<Matrix, 4> drilling_constraint_hessian{
      Matrix{0U, 0U}, Matrix{0U, 0U}, Matrix{0U, 0U}, Matrix{0U, 0U}};
  std::array<Vector, 4> drilling_constraint_third_directional_gradient{
      Vector{0U}, Vector{0U}, Vector{0U}, Vector{0U}};
};

/// @brief Travels with one complete assembly, never accepted frame/history.
struct TrialConstraintData {
  std::filesystem::path source_path;
  double time{0.0};
  std::uint32_t frame_index{0U};
  std::uint32_t increment_index{0U};
  std::uint32_t iteration_count{0U};
  std::vector<DirectionalShellConstraintData> shells;
};

/// @brief Stores a complete trial with independently assembled operators.
struct DynamicAssembly {
  DynamicAcceptedFrameState frame;
  SparseMatrix mass;
  SparseMatrix position_tangent;
  SparseMatrix velocity_tangent;
  SparseMatrix internal_tangent;
  SparseMatrix effective_tangent;
  SparseMatrix drilling_metric;
  Vector force_directional_second;
  Vector force_directional_second_scale;
  std::vector<DynamicRotationKinematics> rotations;
  double momentum_rate_scale;
  double angular_momentum_rate_scale;
  Vector hidden_second{0U};
  Vector hidden_second_scale{0U};
  // Canonical scalar drilling-constraint derivative, indexed by global node.
  std::optional<SparseMatrix> hidden_gradient;
  std::optional<TrialConstraintData> trial_constraints;
};

/// @brief Builds checked shell evaluators in stable source order.
Result<std::vector<RuntimeShell>> BuildRuntimeShells(const Domain& domain,
                                                     const DofManager& dofs);

/// @brief Returns the maximum initial active-element edge length.
Result<double> ReferenceLength(const Domain& domain);

/// @brief Assembles dead loads and their uncancelled absolute row scales.
Result<std::array<Vector, 2>> AssembleLoads(const Domain& domain,
                                            const DofManager& dofs);

/// @brief Evaluates current kinetic and internal quantities without commit.
Result<DynamicAssembly> EvaluateAssembly(
    const Domain& domain, const DofManager& dofs,
    const std::vector<RuntimeShell>& shells, DynamicAcceptedFrameState frame,
    const Vector& load_scale, double time_increment,
    const HhtParameters& parameters, const ParallelFor& parallel_for,
    bool directional_third);

/// @brief Evaluates an increment's initial iterate without committing history.
/// @param frame New-time trial containing the last accepted coordinates.
/// @param predictor Frozen HHT predictor used for both possible starting points.
/// @param old_acceleration Last accepted physical coordinate acceleration.
/// @return The valid initial assembly with smaller free ||D^T E||_infinity,
/// preferring the predictor on a tie. Invalid predictor chart/geometry selects
/// accepted coordinates; rate, mass and arithmetic failures remain failures.
Result<DynamicAssembly> EvaluateNewtonStart(
    const Domain& domain, const DofManager& dofs,
    const std::vector<RuntimeShell>& shells, DynamicAcceptedFrameState frame,
    const HhtPredictor& predictor, const Vector& old_acceleration,
    const Vector& load_scale, double time_increment,
    const HhtParameters& parameters, const ParallelFor& parallel_for);

/// @brief Builds deterministic CSR structure from an ordered local list.
Result<SparseMatrix> MakeSparse(std::size_t rows, std::size_t columns,
                                std::vector<CooContribution> contributions);

/// @brief Builds the scaled free-space operator D^T K D.
Result<SparseMatrix> ScaledFree(const SparseMatrix& matrix,
                                const DofManager& dofs, double length);

/// @brief Checks normwise backward error in the supplied equation space.
Status CheckLinearSolve(const SparseMatrix& matrix, const Vector& rhs,
                        const Vector& solution);

/// @brief Streams checked canonical solves into the caller's infinity norm.
/// @note The caller owns factorization and the initial norm. Each column is
/// discarded after ascending absolute row-sum accumulation; failures preserve norm.
Status AccumulateCheckedCanonicalNorm(const SparseMatrix& matrix,
                                      TangentSolver& solver, double& norm);

/// @brief Finds analytic admissible node null directions without rank cuts.
std::vector<std::size_t> NullNodes(
    const DofManager& dofs,
    const std::vector<DynamicRotationKinematics>& rotations);

/// @brief Contracts a nodal metric with original left-associated products.
/// @note Nonfinite products remain for the caller's existing validity gate.
double OrderedNullMetric(const SparseMatrix& metric, std::size_t node,
                         const Vector3& direction);

/// @brief Solves the positive-mass and differentiated-null initial equations.
///
/// The returned candidate is finite and exact-zero at constrained DOFs. Its
/// complete free-equation N2-R certificate requires reassembly with that
/// acceleration and the caller's source-uncancelled residual scale.
Result<Vector> InitializeAcceleration(const Domain& domain,
                                      const DofManager& dofs,
                                      const DynamicAssembly& initial,
                                      double length,
                                      LinearSolver& positive_solver,
                                      TangentSolver& tangent_solver);

/// @brief Checks mass/regularity and recovers all exact hidden defects.
Status CompleteConstraintEvidence(
    const Domain& domain, const DofManager& dofs,
    const DynamicAssembly& assembly, double length,
    const std::vector<std::size_t>& initial_null_nodes,
    LinearSolver& positive_solver, TangentSolver& tangent_solver,
    DynamicAcceptedFrameState* frame);

/// @brief Recovers physical support balance using volume momentum rates.
Status CompleteEquilibrium(const Domain& domain, const DofManager& dofs,
                           const DynamicAssembly& assembly,
                           DynamicAcceptedFrameState* frame);

}  // namespace hht_internal
}  // namespace fesa
#endif  // FESA_ANALYSIS_HHT_DYNAMIC_INTERNAL_H_