FESADev/src/fesa/analysis/analysis_model.cpp

#include "fesa/analysis/analysis_model.h"

#include <algorithm>
#include <cstddef>
#include <stdexcept>
#include <string>
#include <vector>

namespace fesa {

Result<AnalysisModel> AnalysisModel::Create(const Domain& domain) {
  if (domain.ThermalStep() != nullptr) {
    if (!domain.Steps().Empty()) {
      return Result<AnalysisModel>::Failure(Status::Failure(
          FailureCategory::kModel,
          {{Severity::kError, "mixed-thermal-procedure",
            domain.ThermalStep()->location, "STEP", domain.ThermalStep()->name,
            "Thermal and structural steps cannot coexist."}}));
    }
    return Result<AnalysisModel>::Success(AnalysisModel{domain});
  }
  if (domain.Steps().Empty()) {
    return Result<AnalysisModel>::Failure(
        Status::Failure(FailureCategory::kInput,
                        {{Severity::kError,
                          "invalid-model-cardinality",
                          {domain.SourcePath(), 0U},
                          "STEP",
                          "0",
                          "AnalysisModel requires exactly one static step."}}));
  }
  if (domain.Steps().Size() > 1U) {
    const auto& second_step = domain.Steps()[1];
    return Result<AnalysisModel>::Failure(
        Status::Failure(FailureCategory::kInput,
                        {{Severity::kError, "unsupported-multiple-step",
                          second_step.Location(), "STEP", second_step.Name(),
                          "AnalysisModel does not support multiple steps."}}));
  }
  return Result<AnalysisModel>::Success(AnalysisModel{domain});
}

const Domain& AnalysisModel::GetDomain() const noexcept { return *domain_; }

const StepDefinition& AnalysisModel::Step() const {
  if (is_thermal_) {
    throw std::logic_error{"Thermal AnalysisModel requires ThermalStep()."};
  }
  return domain_->Steps()[0U];
}

bool AnalysisModel::IsThermal() const noexcept { return is_thermal_; }

const ThermalStepDefinition* AnalysisModel::ThermalStep() const noexcept {
  return is_thermal_ ? domain_->ThermalStep() : nullptr;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveElements() const noexcept {
  return active_elements_;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveBeamElements()
    const noexcept {
  return active_beam_elements_;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveMaterials()
    const noexcept {
  return active_materials_;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveProperties()
    const noexcept {
  return active_properties_;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveSections() const noexcept {
  return active_sections_;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveBoundaryConditions()
    const noexcept {
  return active_boundary_conditions_;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveLoads() const noexcept {
  return active_loads_;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveThermalElements()
    const noexcept {
  return active_thermal_elements_;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveThermalNodes()
    const noexcept {
  return active_thermal_nodes_;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveThermalMaterials()
    const noexcept {
  return active_thermal_materials_;
}

const std::vector<EntityIndex>& AnalysisModel::ActiveThermalSections()
    const noexcept {
  return active_thermal_sections_;
}

const std::vector<SurfacePressureLoadDefinition>&
AnalysisModel::ActiveSurfacePressureLoads() const noexcept {
  static const std::vector<SurfacePressureLoadDefinition> kNoSurfacePressures;
  if (is_thermal_) return kNoSurfacePressures;
  return Step().SurfacePressureLoads();
}

AnalysisModel::AnalysisModel(const Domain& domain) : domain_{&domain} {
  if (domain.ThermalStep() != nullptr) {
    is_thermal_ = true;
    std::vector<bool> active_nodes(domain.Nodes().size(), false);
    const std::size_t first_element = domain.Elements().Size() -
                                      domain.ThermalElements().Size();
    for (std::size_t index = 0U; index < domain.ThermalElements().Size();
         ++index) {
      active_thermal_elements_.push_back(
          static_cast<EntityIndex>(first_element + index));
      active_thermal_sections_.push_back(
          domain.ThermalElements()[index].ThermalSectionIndex());
      active_thermal_materials_.push_back(
          domain.ThermalSections()[domain.ThermalElements()[index]
                                       .ThermalSectionIndex()]
              .material_index);
      for (const EntityIndex node :
           domain.ThermalElements()[index].NodeIndices()) {
        active_nodes[node] = true;
      }
    }
    for (std::size_t index = 0U; index < active_nodes.size(); ++index) {
      if (active_nodes[index]) {
        active_thermal_nodes_.push_back(static_cast<EntityIndex>(index));
      }
    }
    std::sort(active_thermal_sections_.begin(), active_thermal_sections_.end());
    active_thermal_sections_.erase(
        std::unique(active_thermal_sections_.begin(),
                    active_thermal_sections_.end()),
        active_thermal_sections_.end());
    std::sort(active_thermal_materials_.begin(),
              active_thermal_materials_.end());
    active_thermal_materials_.erase(
        std::unique(active_thermal_materials_.begin(),
                    active_thermal_materials_.end()),
        active_thermal_materials_.end());
    return;
  }
  std::vector<bool> reachable_materials(domain.Materials().Size(), false);
  std::vector<bool> reachable_properties(domain.Properties().Size(), false);
  std::vector<bool> reachable_sections(domain.Sections().Size(), false);

  for (std::size_t index = 0U; index < domain.Elements().Size(); ++index) {
    const auto& element = domain.Elements()[index];
    active_elements_.push_back(static_cast<EntityIndex>(index));
    reachable_materials[element.MaterialIndex()] = true;
    reachable_properties[element.PropertyIndex()] = true;
  }

  for (std::size_t index = 0U; index < domain.BeamElements().Size(); ++index) {
    const auto& element = domain.BeamElements()[index];
    active_beam_elements_.push_back(static_cast<EntityIndex>(index));
    reachable_sections[element.section_index] = true;
  }

  // Ascending vector positions are the stable internal order, independent
  // of first reachability or duplicate element assignments.
  for (std::size_t index = 0U; index < reachable_materials.size(); ++index) {
    if (reachable_materials[index]) {
      active_materials_.push_back(static_cast<EntityIndex>(index));
    }
  }
  for (std::size_t index = 0U; index < reachable_properties.size(); ++index) {
    if (reachable_properties[index]) {
      active_properties_.push_back(static_cast<EntityIndex>(index));
    }
  }
  for (std::size_t index = 0U; index < reachable_sections.size(); ++index) {
    if (reachable_sections[index]) {
      active_sections_.push_back(static_cast<EntityIndex>(index));
    }
  }

  for (std::size_t index = 0U; index < Step().BoundaryConditions().size();
       ++index) {
    active_boundary_conditions_.push_back(static_cast<EntityIndex>(index));
  }
  for (std::size_t index = 0U; index < Step().Loads().size(); ++index) {
    active_loads_.push_back(static_cast<EntityIndex>(index));
  }
}

}  // namespace fesa