FESADev/src/fesa/analysis/analysis_state.cpp

#include "fesa/analysis/analysis_state.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <limits>
#include <string>
#include <utility>

namespace fesa {
namespace {

constexpr std::size_t kShellLocationsPerElement = 4U;

/// @brief Creates a structured failure without mutating existing state.
Status ShellCandidateFailure(const std::string& code,
                             const std::string& identity,
                             const std::string& message) {
  return Status::Failure(
      FailureCategory::kModel,
      {{Severity::kError, code, {}, "ANALYSIS_STATE", identity, message}});
}

template <std::size_t Size>
/// @brief Tests one fixed-size candidate component array for finite values.
bool IsFinite(const std::array<double, Size>& values) {
  return std::all_of(values.begin(), values.end(),
                     [](const double value) { return std::isfinite(value); });
}

/// @brief Tests every physical shell row component for finite values.
bool IsFinite(const ShellResultRow& row) {
  if (!IsFinite(row.natural_coordinates) || !IsFinite(row.generalized_strain) ||
      !IsFinite(row.section_resultant)) {
    return false;
  }
  for (const auto& axis : row.local_frame) {
    if (!IsFinite(axis)) {
      return false;
    }
  }
  return std::all_of(row.stress.begin(), row.stress.end(),
                     [](const ShellSectionStressRow& stress) {
                       return std::isfinite(stress.zeta) &&
                              IsFinite(stress.components);
                     });
}

}  // namespace

AnalysisState AnalysisState::Create(const DofManager& dofs,
                                    StepFrameIdentity identity) {
  return AnalysisState{dofs.FullDofCount(), std::move(identity)};
}

Vector& AnalysisState::Displacement() noexcept { return displacement_; }

const Vector& AnalysisState::Displacement() const noexcept {
  return displacement_;
}

Vector& AnalysisState::ExternalForce() noexcept { return external_force_; }

const Vector& AnalysisState::ExternalForce() const noexcept {
  return external_force_;
}

Vector& AnalysisState::InternalForce() noexcept { return internal_force_; }

const Vector& AnalysisState::InternalForce() const noexcept {
  return internal_force_;
}

Vector& AnalysisState::Residual() noexcept { return residual_; }

const Vector& AnalysisState::Residual() const noexcept { return residual_; }

Vector& AnalysisState::Reaction() noexcept { return reaction_; }

const Vector& AnalysisState::Reaction() const noexcept { return reaction_; }

const StepFrameIdentity& AnalysisState::Identity() const noexcept {
  return identity_;
}

std::vector<EndpointResultRow>& AnalysisState::EndpointResults() noexcept {
  return endpoint_results_;
}

const std::vector<EndpointResultRow>& AnalysisState::EndpointResults()
    const noexcept {
  return endpoint_results_;
}

std::vector<GaussResultRow>& AnalysisState::GaussResults() noexcept {
  return gauss_results_;
}

const std::vector<GaussResultRow>& AnalysisState::GaussResults()
    const noexcept {
  return gauss_results_;
}

std::vector<StressS11Row>& AnalysisState::StressResults() noexcept {
  return stress_results_;
}

const std::vector<StressS11Row>& AnalysisState::StressResults() const noexcept {
  return stress_results_;
}

Status AnalysisState::CommitShellResults(
    const std::vector<EntityIndex>& expected_element_order,
    ShellStateCandidate candidate) {
  if (expected_element_order.size() >
      (std::numeric_limits<std::size_t>::max)() / kShellLocationsPerElement) {
    return ShellCandidateFailure(
        "invalid-shell-state-inventory", identity_.step_name,
        "The expected shell result inventory is too large.");
  }
  const std::size_t expected_row_count =
      expected_element_order.size() * kShellLocationsPerElement;
  if (candidate.rows.size() != expected_row_count) {
    return ShellCandidateFailure(
        "invalid-shell-state-inventory", identity_.step_name,
        "Shell results require exactly four rows per expected element.");
  }
  if (std::adjacent_find(expected_element_order.begin(),
                         expected_element_order.end(),
                         [](const EntityIndex left, const EntityIndex right) {
                           return left >= right;
                         }) != expected_element_order.end()) {
    return ShellCandidateFailure(
        "invalid-shell-state-inventory", identity_.step_name,
        "Expected shell elements must be unique and in stable index order.");
  }

  const std::array<ShellMidsurfaceLocation, kShellLocationsPerElement>
      expected_locations{
          ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2,
          ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4};
  const double gauss = 1.0 / std::sqrt(3.0);
  const std::array<std::array<double, 2>, kShellLocationsPerElement>
      expected_coordinates{std::array<double, 2>{-gauss, -gauss},
                           std::array<double, 2>{gauss, -gauss},
                           std::array<double, 2>{gauss, gauss},
                           std::array<double, 2>{-gauss, gauss}};
  const std::array<ShellSectionPosition, 3> expected_positions{
      ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle,
      ShellSectionPosition::kTop};
  constexpr std::array<double, 3> expected_zeta{-1.0, 0.0, 1.0};
  for (std::size_t element_order = 0U;
       element_order < expected_element_order.size(); ++element_order) {
    for (std::size_t point = 0U; point < kShellLocationsPerElement; ++point) {
      const auto& row =
          candidate.rows[element_order * kShellLocationsPerElement + point];
      if (row.element != expected_element_order[element_order] ||
          row.location != expected_locations[point] ||
          row.natural_coordinates != expected_coordinates[point]) {
        return ShellCandidateFailure(
            "invalid-shell-state-inventory", std::to_string(row.element),
            "Shell rows must preserve element and GP1 through GP4 identity.");
      }
      for (std::size_t position = 0U; position < expected_positions.size();
           ++position) {
        if (row.stress[position].position != expected_positions[position] ||
            row.stress[position].zeta != expected_zeta[position]) {
          return ShellCandidateFailure(
              "invalid-shell-state-inventory", std::to_string(row.element),
              "Shell stress rows require BOTTOM, MIDDLE, TOP identity.");
        }
      }
      if (!IsFinite(row)) {
        return ShellCandidateFailure(
            "nonfinite-shell-state-value", std::to_string(row.element),
            "Shell result rows must contain only finite values.");
      }
    }
  }

  if (!std::isfinite(candidate.physical_strain_energy) ||
      !IsFinite(candidate.equilibrium) ||
      !IsFinite(candidate.verification_metrics)) {
    return ShellCandidateFailure(
        "nonfinite-shell-state-value", identity_.step_name,
        "Shell energy, equilibrium, and normalized metrics must be finite.");
  }

  shell_results_ = std::move(candidate.rows);
  physical_strain_energy_ = candidate.physical_strain_energy;
  equilibrium_ = candidate.equilibrium;
  verification_metrics_ = candidate.verification_metrics;
  return Status::Ok();
}

const std::vector<ShellResultRow>& AnalysisState::ShellResults()
    const noexcept {
  return shell_results_;
}

double AnalysisState::PhysicalStrainEnergy() const noexcept {
  return physical_strain_energy_;
}

const std::array<double, 6>& AnalysisState::Equilibrium() const noexcept {
  return equilibrium_;
}

const std::array<double, 3>& AnalysisState::VerificationMetrics()
    const noexcept {
  return verification_metrics_;
}

AnalysisState::AnalysisState(std::size_t full_dof_count,
                             StepFrameIdentity identity)
    : identity_{std::move(identity)},
      displacement_{full_dof_count},
      external_force_{full_dof_count},
      internal_force_{full_dof_count},
      residual_{full_dof_count},
      reaction_{full_dof_count} {}

}  // namespace fesa