#include <gtest/gtest.h>

#include <algorithm>
#include <filesystem>
#include <fstream>
#include <functional>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>

#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/fem/dof_manager.h"
#include "fesa/io/abaqus/domain_mapper.h"
#include "fesa/io/abaqus/input_reader.h"

namespace {
std::string Deck(const std::string& procedure = "*Static\n1,1,1,1\n",
                 const std::string& boundary = "10,1,3\n",
                 const std::string& extra = "") {
  return "*Part,name=P\n*Node\n10,0,0,0\n20,1,0,0\n30,1,1,0\n40,0,1,0\n"
         "50,0,0,1\n60,1,0,1\n70,1,1,1\n80,0,1,1\n"
         "*Element,type=C3D8,elset=E\n99,10,20,30,40,50,60,70,80\n"
         "*Solid Section,elset=E,material=M\n,\n*End Part\n"
         "*Assembly,name=A\n*Instance,name=I,part=P\n*End Instance\n"
         "*End Assembly\n*Material,name=M\n*Elastic\n1000,.3\n" +
         extra + "*Step,name=Solid,inc=100\n" + procedure +
         (boundary.empty() ? "" : "*Boundary\n" + boundary) + "*End Step\n";
}
fesa::Result<fesa::Domain> Map(const std::string& text) {
  const auto path =
      std::filesystem::temp_directory_path() / "fesa-hex8-semantics.inp";
  {
    std::ofstream out{path};
    out << text;
  }
  auto syntax = fesa::AbaqusInputReader{}.Read(path);
  std::filesystem::remove(path);
  if (!syntax.HasValue())
    return fesa::Result<fesa::Domain>::Failure(syntax.GetStatus());
  return fesa::AbaqusDomainMapper{}.Map(syntax.Value());
}
TEST(Hex8Semantics, MapsEightSourceNodesAndExactlyThreePhysicalDofs) {
  auto domain = Map(Deck());
  ASSERT_TRUE(domain.HasValue())
      << (domain.GetStatus().Diagnostics().empty()
              ? ""
              : domain.GetStatus().Diagnostics()[0].message);
  ASSERT_EQ(domain.Value().Elements().Size(), 1U);
  EXPECT_EQ(domain.Value().Elements()[0].SourceElementType(), "C3D8");
  EXPECT_EQ(domain.Value().Elements()[0].SourceId().source_label, 99);
  EXPECT_EQ(domain.Value().Elements()[0].NodeIndices().size(), 8U);
  auto model = fesa::AnalysisModel::Create(domain.Value());
  ASSERT_TRUE(model.HasValue());
  auto dofs = fesa::DofManager::Create(model.Value());
  ASSERT_TRUE(dofs.HasValue());
  EXPECT_EQ(dofs.Value().FullDofCount(), 24U);
  EXPECT_EQ(dofs.Value().ConstrainedDofCount(), 3U);
}
TEST(Hex8Semantics, InactiveZeroRotationsWarnAndDoNotCreateConstraints) {
  auto domain = Map(Deck("*Static\n1,1,1,1\n", "10,1,6\n"));
  ASSERT_TRUE(domain.HasValue())
      << (domain.GetStatus().Diagnostics().empty()
              ? ""
              : domain.GetStatus().Diagnostics()[0].message);
  EXPECT_EQ(domain.Value().Steps()[0].BoundaryConditions().size(), 3U);
  ASSERT_FALSE(domain.Value().Warnings().empty());
  EXPECT_EQ(domain.Value().Warnings().back().keyword, "BOUNDARY");
  EXPECT_EQ(domain.Value().Warnings().back().entity_identity, "10");
  EXPECT_GT(domain.Value().Warnings().back().location.line, 0U);
}
TEST(Hex8Semantics, RejectsNonzeroRotationsAndAllRotationalLoads) {
  EXPECT_FALSE(Map(Deck("*Static\n1,1,1,1\n", "10,4,4,1\n")).HasValue());
  EXPECT_FALSE(
      Map(Deck("*Static\n1,1,1,1\n", "10,1,3\n*Cload\n20,4,0\n")).HasValue());
}

std::string Replace(std::string input, const std::string& from,
                    const std::string& to) {
  const auto position = input.find(from);
  if (position == std::string::npos)
    throw std::logic_error{"Missing mutation source"};
  input.replace(position, from.size(), to);
  return input;
}
TEST(Hex8Semantics, RotationNoOpPreservesPhysicalMapAndLoads) {
  auto physical =
      Map(Deck("*Static\n1,1,1,1\n", "10,1,3,0\n*Cload\n20,2,17\n"));
  auto extra =
      Map(Deck("*Static\n1,1,1,1\n", "10,1,6,0\n20,4,6,-0\n*Cload\n20,2,17\n"));
  ASSERT_TRUE(physical.HasValue());
  ASSERT_TRUE(extra.HasValue());
  auto pm = fesa::AnalysisModel::Create(physical.Value());
  auto em = fesa::AnalysisModel::Create(extra.Value());
  auto pd = fesa::DofManager::Create(pm.Value());
  auto ed = fesa::DofManager::Create(em.Value());
  ASSERT_TRUE(pd.HasValue());
  ASSERT_TRUE(ed.HasValue());
  EXPECT_EQ(pd.Value().FreeDofs(), ed.Value().FreeDofs());
  EXPECT_EQ(pd.Value().ConstrainedDofs(), ed.Value().ConstrainedDofs());
  EXPECT_EQ(pd.Value().GetSparsePattern().column_indices,
            ed.Value().GetSparsePattern().column_indices);
  EXPECT_EQ(ed.Value().Profile(), fesa::DofProfile::kSolidStructural);
  EXPECT_EQ(ed.Value().ComponentsPerNode(), 3U);
  EXPECT_THROW(ed.Value().FullDof(0, fesa::DofComponent::kUrx),
               std::out_of_range);
  const auto& element = extra.Value().Elements()[0];
  auto scatter = ed.Value().ElementScatter(
      {element.SourceId(),
       element.NodeIndices(),
       {fesa::DofComponent::kUx, fesa::DofComponent::kUy,
        fesa::DofComponent::kUz}});
  ASSERT_TRUE(scatter.HasValue());
  EXPECT_EQ(scatter.Value().size(), 24U);
  EXPECT_EQ(scatter.Value().front(), 0U);
  EXPECT_EQ(scatter.Value().back(), 23U);
  auto pf = fesa::LoadAssembler::AssembleFullNodalLoad(pm.Value(), pd.Value());
  auto ef = fesa::LoadAssembler::AssembleFullNodalLoad(em.Value(), ed.Value());
  ASSERT_TRUE(pf.HasValue());
  ASSERT_TRUE(ef.HasValue());
  for (std::size_t i = 0; i < 24U; ++i) EXPECT_EQ(pf.Value()[i], ef.Value()[i]);
}
TEST(Hex8Semantics, RejectsTinyRotationBadTargetsOptionsAndDofRanges) {
  for (const std::string row :
       {"10,4,6,1e-100\n", "Missing,4,6,0\n", "10,0,3\n", "10,1,7\n"})
    EXPECT_FALSE(Map(Deck("*Static\n1,1,1,1\n", row)).HasValue()) << row;
  for (const std::string option : {"type=VELOCITY", "amplitude=A", "op=NEW"})
    EXPECT_FALSE(
        Map(Replace(Deck(), "*Boundary", "*Boundary," + option)).HasValue());
  for (const std::string row : {"10,4,6\n", "10,4,6,0\n", "10,4,6,\n"}) {
    auto result = Map(Deck("*Static\n1,1,1,1\n", row));
    ASSERT_TRUE(result.HasValue()) << row;
    EXPECT_TRUE(result.Value().Steps()[0].BoundaryConditions().empty());
    ASSERT_EQ(result.Value().Warnings().size(), 1U);
    EXPECT_EQ(result.Value().Warnings()[0].code,
              "inactive-solid-rotation-boundary");
    EXPECT_EQ(result.Value().Warnings()[0].severity, fesa::Severity::kWarning);
  }
}
TEST(Hex8Semantics, AllowsEmptyFreeAndNoBoundaryModels) {
  std::string all;
  for (int node = 10; node <= 80; node += 10)
    all += std::to_string(node) + ",1,6\n";
  for (const auto& boundaries : {all, std::string{}}) {
    auto domain = Map(Deck("*Static\n1,1,1,1\n", boundaries));
    ASSERT_TRUE(domain.HasValue());
    auto model = fesa::AnalysisModel::Create(domain.Value());
    auto dofs = fesa::DofManager::Create(model.Value());
    ASSERT_TRUE(dofs.HasValue());
    EXPECT_EQ(dofs.Value().FullDofCount(), 24U);
    EXPECT_EQ(dofs.Value().FreeDofCount(), boundaries.empty() ? 24U : 0U);
  }
}
TEST(Hex8Semantics, RetainsStaticDensityTargetsAndDirectSchedule) {
  auto domain =
      Map(Deck("*Static,direct\n.1,1,\n", "10,1,3,.25\n", "*Density\n2\n"));
  ASSERT_TRUE(domain.HasValue())
      << (domain.GetStatus().Diagnostics().empty()
              ? ""
              : domain.GetStatus().Diagnostics()[0].message);
  EXPECT_TRUE(domain.Value().IsSolid());
  ASSERT_EQ(domain.Value().SolidSections().Size(), 1U);
  EXPECT_EQ(domain.Value().SolidSections()[0].MaterialIndex(), 0U);
  EXPECT_EQ(domain.Value().Densities()[0].rho, 2.0);
  const auto& step = domain.Value().Steps()[0];
  EXPECT_EQ(step.Name(), "Solid");
  EXPECT_TRUE(step.IsDirectStatic());
  EXPECT_EQ(step.IncrementCount(), 10U);
  EXPECT_EQ(step.MaximumIncrements(), 100U);
  auto model = fesa::AnalysisModel::Create(domain.Value());
  auto dofs = fesa::DofManager::Create(model.Value());
  EXPECT_EQ(dofs.Value().PrescribedValues()[0], .25);
  for (const std::string values : {".3,1\n", "0,1\n", ".1,11\n"})
    EXPECT_FALSE(Map(Deck("*Static,direct\n" + values)).HasValue());
  EXPECT_FALSE(Map(Replace(Deck(), "1000,.3", "1000,.5")).HasValue());
  EXPECT_FALSE(Map(Replace(Deck(), "1000,.3", "1000,-1")).HasValue());
}
TEST(Hex8Semantics, SamplesIndependentAmplitudeLoadsAtActualEndpoints) {
  const std::string amplitudes =
      "*Density\n2\n*Amplitude,name=A\n0,0,.1,1\n10,1\n"
      "*Amplitude,name=B,definition=TABULAR,time=STEP "
      "TIME,value=RELATIVE\n0,2,1,4\n";
  const auto text = Deck("*Dynamic,direct\n.1,10,\n",
                         "10,1,6\n*Cload,amplitude=a\n20,1,10\n"
                         "*Cload,amplitude=B\n20,1,3\n*Cload\n20,1,5\n",
                         amplitudes);
  auto domain = Map(text);
  ASSERT_TRUE(domain.HasValue())
      << (domain.GetStatus().Diagnostics().empty()
              ? ""
              : domain.GetStatus().Diagnostics()[0].message);
  ASSERT_EQ(domain.Value().Amplitudes().size(), 2U);
  const auto& amplitude = domain.Value().Amplitudes()[0];
  EXPECT_EQ(amplitude.name_raw, "A");
  EXPECT_EQ(amplitude.Evaluate(0), 0);
  EXPECT_EQ(amplitude.Evaluate(.05), .5);
  EXPECT_EQ(amplitude.Evaluate(.1), 1);
  EXPECT_EQ(amplitude.Evaluate(20), 1);
  EXPECT_FALSE(domain.Value().Dynamic()->nlgeom);
  EXPECT_EQ(domain.Value().Dynamic()->increment_count, 100U);
  const auto& loads = domain.Value().Steps()[0].ConcentratedLoads();
  EXPECT_EQ(loads[0].AmplitudeIndex(), std::optional<fesa::EntityIndex>{0U});
  EXPECT_EQ(loads[1].AmplitudeIndex(), std::optional<fesa::EntityIndex>{1U});
  EXPECT_EQ(loads[0].AmplitudeName(), "a");
  EXPECT_GT(loads[0].AmplitudeLocation().line, 0U);
  auto model = fesa::AnalysisModel::Create(domain.Value());
  auto dofs = fesa::DofManager::Create(model.Value());
  ASSERT_TRUE(dofs.HasValue());
  for (double time : {0., .05, .1, 1., 10.}) {
    auto load = fesa::LoadAssembler::AssembleFullNodalLoad(model.Value(),
                                                           dofs.Value(), time);
    ASSERT_TRUE(load.HasValue());
    EXPECT_DOUBLE_EQ(load.Value()[3],
                     amplitude.Evaluate(time) * 10 +
                         domain.Value().Amplitudes()[1].Evaluate(time) * 3 + 5);
  }
  EXPECT_FALSE(Map(Replace(text, "10,1,6", "10,1,3,1e-100")).HasValue());
  EXPECT_FALSE(
      Map(Replace(text, "amplitude=a", "amplitude=missing")).HasValue());
  EXPECT_TRUE(
      Map(Replace(text, "name=Solid,inc=100", "name=Solid,nlgeom=YES,inc=100"))
          .HasValue());
}
TEST(Hex8Semantics, RejectsAmplitudeOptionsDataAndCrossFamilyUse) {
  const std::string prefix = "*Density\n2\n";
  for (const std::string amplitude :
       {"*Amplitude,name=A,time=TOTAL TIME\n0,0\n",
        "*Amplitude,name=A,definition=SMOOTH STEP\n0,0\n",
        "*Amplitude,name=A,value=ABSOLUTE\n0,0\n",
        "*Amplitude,name=A,smooth=0\n0,0\n", "*Amplitude,name=A\n0,0,1\n",
        "*Amplitude,name=A\n0,0,0,1\n", "*Amplitude,name=A\n1,0,0,1\n",
        "*Amplitude,name=A\n-1,0\n", "*Amplitude,name=A\n0,nan\n",
        "*Amplitude,name=A\n",
        "*Amplitude,name=A\n0,0\n*Amplitude,name=a\n0,1\n"})
    EXPECT_FALSE(
        Map(Deck("*Dynamic,direct\n.1,1\n", "10,1,3\n", prefix + amplitude))
            .HasValue())
        << amplitude;
  EXPECT_FALSE(
      Map(Deck("*Static\n1,1,1,1\n", "10,1,3\n", "*Amplitude,name=A\n0,0\n"))
          .HasValue());
  EXPECT_FALSE(Map(Deck("*Dynamic,direct\n.1,1\n", "10,1,3\n")).HasValue());
}
TEST(Hex8Semantics, RejectsAliasesConnectivityThicknessAndMixedModels) {
  for (const std::string type : {"C3D8R", "C3D8I", "C3D8H", "C3D20"})
    EXPECT_FALSE(Map(Replace(Deck(), "type=C3D8", "type=" + type)).HasValue());
  EXPECT_FALSE(Map(Replace(Deck(), "99,10,20,30,40,50,60,70,80",
                           "99,10,20,30,40,50,60,70,70"))
                   .HasValue());
  EXPECT_FALSE(Map(Replace(Deck(), "99,10,20,30,40,50,60,70,80",
                           "99,10,20,30,40,50,60,70,999"))
                   .HasValue());
  EXPECT_FALSE(
      Map(Replace(Deck(), "material=M\n,", "material=M\n2")).HasValue());
  EXPECT_FALSE(Map(Replace(Deck(), "*Solid Section",
                           "*Element,type=S4\n100,10,20,30,40\n*Solid Section"))
                   .HasValue());
}
TEST(Hex8Semantics, PreservesShellSixComponentProfile) {
  auto text = Replace(Deck(), "50,0,0,1\n60,1,0,1\n70,1,1,1\n80,0,1,1\n", "");
  text = Replace(text, "*Element,type=C3D8,elset=E\n99,10,20,30,40,50,60,70,80",
                 "*Element,type=S4\n99,10,20,30,40\n*Elset,elset=E\n99");
  text = Replace(text, "*Solid Section,elset=E,material=M\n,",
                 "*Shell Section,elset=E,material=M\n.1");
  auto domain = Map(text);
  ASSERT_TRUE(domain.HasValue());
  auto model = fesa::AnalysisModel::Create(domain.Value());
  auto dofs = fesa::DofManager::Create(model.Value());
  ASSERT_TRUE(dofs.HasValue());
  EXPECT_EQ(dofs.Value().Profile(), fesa::DofProfile::kStructural);
  EXPECT_EQ(dofs.Value().ComponentsPerNode(), 6U);
  EXPECT_EQ(dofs.Value().FullDofCount(), 24U);
  EXPECT_EQ(dofs.Value().FullDof(1, fesa::DofComponent::kUrx), 9U);
}
TEST(Hex8Semantics, MapsAllFiveUnmodifiedRequiredSourceDecks) {
  auto root = std::filesystem::path{__FILE__}.parent_path();
  for (int i = 0; i < 4; ++i) root = root.parent_path();
  for (const std::string name :
       {"hexa_linearstatic.inp", "hexa_nonlinearstatic.inp", "hexa_freq.inp",
        "hexa_dynlinear.inp", "hexa_dynnonlinear.inp"}) {
    auto syntax =
        fesa::AbaqusInputReader{}.Read(root / "reference/hexa" / name);
    ASSERT_TRUE(syntax.HasValue()) << name;
    auto domain = fesa::AbaqusDomainMapper{}.Map(syntax.Value());
    ASSERT_TRUE(domain.HasValue())
        << name << ": "
        << (domain.GetStatus().Diagnostics().empty()
                ? ""
                : domain.GetStatus().Diagnostics()[0].message);
    EXPECT_EQ(domain.Value().SolidElements().Size(), 300U);
    EXPECT_EQ(domain.Value().Nodes().size(), 442U);
    EXPECT_EQ(std::count_if(domain.Value().Warnings().begin(),
                            domain.Value().Warnings().end(),
                            [](const auto& warning) {
                              return warning.code ==
                                     "inactive-solid-rotation-boundary";
                            }),
              3);
    auto model = fesa::AnalysisModel::Create(domain.Value());
    auto dofs = fesa::DofManager::Create(model.Value());
    ASSERT_TRUE(dofs.HasValue()) << name;
    EXPECT_EQ(dofs.Value().FullDofCount(), 1326U);
    EXPECT_EQ(dofs.Value().ConstrainedDofCount(), 51U);
  }
}

class ConstantMatrixElement final : public fesa::Element {
 public:
  ConstantMatrixElement(const fesa::ElementDefinition& definition, double entry)
      : layout_{definition.SourceId(),
                definition.NodeIndices(),
                {fesa::DofComponent::kUx, fesa::DofComponent::kUy,
                 fesa::DofComponent::kUz}},
        entry_{entry} {}
  const fesa::ElementDofLayout& DofLayout() const noexcept override {
    return layout_;
  }
  fesa::Result<fesa::ElementStiffnessContribution> ComputeStiffness()
      const override {
    fesa::Matrix values{24U, 24U};
    for (std::size_t row = 0U; row < 24U; ++row)
      for (std::size_t column = 0U; column < 24U; ++column)
        values(row, column) = entry_;
    return fesa::Result<fesa::ElementStiffnessContribution>::Success(
        {layout_, std::move(values)});
  }
  fesa::Result<fesa::ElementResultBundle> Recover(
      const fesa::Vector&) const override {
    return fesa::Result<fesa::ElementResultBundle>::Failure(
        fesa::Status::Failure(fesa::FailureCategory::kModel,
                              {{fesa::Severity::kError,
                                "test-no-recovery",
                                {},
                                "",
                                "",
                                "Assembly fixture has no recovery."}}));
  }

 private:
  fesa::ElementDofLayout layout_;
  double entry_;
};
TEST(Hex8Semantics, AssemblesFullSymmetricSolidMatricesInStableSourceOrder) {
  const auto text =
      Replace(Deck(), "99,10,20,30,40,50,60,70,80\n",
              "99,10,20,30,40,50,60,70,80\n100,10,20,30,40,50,60,70,80\n"
              "101,10,20,30,40,50,60,70,80\n");
  auto domain = Map(text);
  ASSERT_TRUE(domain.HasValue());
  auto model = fesa::AnalysisModel::Create(domain.Value());
  ConstantMatrixElement first{domain.Value().Elements()[0], 1e16};
  ConstantMatrixElement second{domain.Value().Elements()[1], -1e16};
  ConstantMatrixElement third{domain.Value().Elements()[2], 1};
  const fesa::ElementView elements{std::cref(first), std::cref(second),
                                   std::cref(third)};
  fesa::DofManager dofs;
  ASSERT_TRUE(dofs.Build(model.Value(), elements).IsOk());
  auto serial = fesa::SparseAssembler::Assemble(elements, dofs,
                                                fesa::SerialParallelFor{});
  auto parallel =
      fesa::SparseAssembler::Assemble(elements, dofs, fesa::TbbParallelFor{});
  ASSERT_TRUE(serial.HasValue());
  ASSERT_TRUE(parallel.HasValue());
  EXPECT_EQ(serial.Value().Rows(), 24U);
  EXPECT_EQ(serial.Value().Columns(), 24U);
  EXPECT_EQ(serial.Value().Values().size(), 576U);
  for (double value : serial.Value().Values()) EXPECT_EQ(value, 1.0);
  EXPECT_EQ(serial.Value().RowOffsets(), parallel.Value().RowOffsets());
  EXPECT_EQ(serial.Value().ColumnIndices(), parallel.Value().ColumnIndices());
  EXPECT_EQ(serial.Value().Values(), parallel.Value().Values());
}
TEST(Hex8Semantics, ValidatesOutputNoOpsAndRejectsUnapprovedDynamicChildren) {
  for (const std::string output :
       {"*Output,field,variable=PRESELECT\n*Node Output\nU,RF\n*Element "
        "Output,directions=YES\nS,E\n",
        "*Output,history,variable=PRESELECT\n"})
    EXPECT_TRUE(
        Map(Replace(Deck(), "*End Step", output + "*End Step")).HasValue());
  EXPECT_FALSE(Map(Replace(Deck(), "*End Step",
                           "*Output,field,variable=OTHER\n*End Step"))
                   .HasValue());
  EXPECT_FALSE(
      Map(Replace(
              Deck(), "*End Step",
              "*Output,field\n*Element Output,position=NODES\nS\n*End Step"))
          .HasValue());
  auto dynamic = Deck("*Dynamic,direct\n.1,1\n", "10,1,3\n", "*Density\n2\n");
  EXPECT_FALSE(
      Map(Replace(
              dynamic, "*End Step",
              "*Output,field,variable=PRESELECT\n*Node Output\nU\n*End Step"))
          .HasValue());
  EXPECT_FALSE(Map(Replace(dynamic, "*End Step",
                           "*Restart,write,frequency=1\n*End Step"))
                   .HasValue());
}

fesa::ModelDefinition SolidDefinition() {
  fesa::ModelDefinition definition;
  const std::array<std::array<double, 3>, 8> coordinates{{{0, 0, 0},
                                                          {1, 0, 0},
                                                          {1, 1, 0},
                                                          {0, 1, 0},
                                                          {0, 0, 1},
                                                          {1, 0, 1},
                                                          {1, 1, 1},
                                                          {0, 1, 1}}};
  for (std::size_t i = 0; i < coordinates.size(); ++i)
    definition.nodes.push_back({{"I", static_cast<std::int64_t>(10 * (i + 1)),
                                 std::to_string(10 * (i + 1))},
                                coordinates[i],
                                {"solid.inp", i + 1}});
  definition.materials.emplace_back("M", 1000, .3,
                                    fesa::SourceLocation{"solid.inp", 20});
  definition.solid_sections.emplace_back(fesa::SourceEntityId{"P", 1, "E"}, "E",
                                         0,
                                         fesa::SourceLocation{"solid.inp", 21});
  definition.solid_elements.emplace_back(
      fesa::SourceEntityId{"I", 99, "99"},
      std::array<fesa::EntityIndex, 8>{{0, 1, 2, 3, 4, 5, 6, 7}}, 0, 0,
      fesa::SourceLocation{"solid.inp", 22});
  definition.steps.emplace_back(
      "Solid", std::vector<fesa::PrescribedDisplacementDefinition>{},
      std::vector<fesa::NodalLoad>{}, 1., 1., 1., 1.,
      fesa::SourceLocation{"solid.inp", 23});
  return definition;
}
TEST(Hex8Semantics, RejectsInvalidProgrammaticSolidCandidates) {
  EXPECT_TRUE(fesa::Domain::Create(SolidDefinition()).HasValue());
  auto repeated = SolidDefinition();
  repeated.solid_elements[0].node_indices[7] = 0;
  EXPECT_FALSE(fesa::Domain::Create(std::move(repeated)).HasValue());
  auto material = SolidDefinition();
  material.materials[0].poisson_ratio = .5;
  EXPECT_FALSE(fesa::Domain::Create(std::move(material)).HasValue());
  auto section = SolidDefinition();
  section.solid_elements[0].section_index = 99;
  EXPECT_FALSE(fesa::Domain::Create(std::move(section)).HasValue());
  auto rotational = SolidDefinition();
  rotational.steps[0].loads.push_back({"10", 4, 0, {"solid.inp", 24}});
  EXPECT_FALSE(fesa::Domain::Create(std::move(rotational)).HasValue());
  auto unused = SolidDefinition();
  unused.nodes.push_back({{"I", 90, "90"}, {2, 2, 2}, {"solid.inp", 25}});
  EXPECT_FALSE(fesa::Domain::Create(std::move(unused)).HasValue());
}

TEST(Hex8Semantics, RejectsSecondProcedureBeforeReadingItsOptionsOrData) {
  const std::array<std::string, 4> procedures{
      "*Static\n1,1,1,1\n", "*Static,direct\n.1,1\n", "*Frequency\n2\n",
      "*Dynamic,direct\n.1,1\n"};
  for (const auto& first : procedures) {
    ASSERT_TRUE(Map(Deck(first, "", "*Density\n2\n")).HasValue()) << first;
    for (const auto& second : procedures) {
      SCOPED_TRACE(first + second);
      const auto result = Map(Deck(first + second, "", "*Density\n2\n"));
      ASSERT_FALSE(result.HasValue());
      ASSERT_FALSE(result.GetStatus().Diagnostics().empty());
      EXPECT_EQ(result.GetStatus().Diagnostics()[0].code,
                "unsupported-analysis-procedure");
    }
    const auto malformed =
        Map(Deck(first + "*Static,direct,unapproved=YES\nnot-a-number\n", "",
                 "*Density\n2\n"));
    ASSERT_FALSE(malformed.HasValue());
    ASSERT_FALSE(malformed.GetStatus().Diagnostics().empty());
    EXPECT_EQ(malformed.GetStatus().Diagnostics()[0].code,
              "unsupported-analysis-procedure");
    EXPECT_EQ(malformed.GetStatus().Diagnostics()[0].keyword, "STATIC");
  }
}
TEST(Hex8Semantics, AllowsOnlyApprovedStaticRestartNoOp) {
  auto approved = Map(
      Replace(Deck(), "*End Step", "*Restart,write,frequency=0\n*End Step"));
  ASSERT_TRUE(approved.HasValue());
  ASSERT_EQ(approved.Value().Warnings().size(), 1U);
  EXPECT_EQ(approved.Value().Warnings()[0].keyword, "RESTART");
  for (const std::string options :
       {"read", "read,step=1,inc=1", "write", "write,frequency=1",
        "write,frequency=0,endstep", "write,frequency=0,numberinterval=1",
        "write=YES,frequency=0"}) {
    SCOPED_TRACE(options);
    auto result = Map(
        Replace(Deck(), "*End Step", "*Restart," + options + "\n*End Step"));
    ASSERT_FALSE(result.HasValue());
    ASSERT_FALSE(result.GetStatus().Diagnostics().empty());
    EXPECT_EQ(result.GetStatus().Diagnostics()[0].keyword, "RESTART");
  }
}
TEST(Hex8Semantics,
     PreservesShuffledSourceLabelsAndRejectsDuplicateIdentities) {
  const auto deck = Replace(Deck(),
                            "10,0,0,0\n20,1,0,0\n30,1,1,0\n40,0,1,0\n"
                            "50,0,0,1\n60,1,0,1\n70,1,1,1\n80,0,1,1\n",
                            "80,0,1,1\n20,1,0,0\n60,1,0,1\n40,0,1,0\n"
                            "10,0,0,0\n70,1,1,1\n30,1,1,0\n50,0,0,1\n");
  auto domain = Map(deck);
  ASSERT_TRUE(domain.HasValue());
  EXPECT_EQ(domain.Value().Nodes()[0].source_id.source_label, 80);
  EXPECT_EQ(domain.Value().Nodes()[4].source_id.source_label, 10);
  EXPECT_EQ(domain.Value().Elements()[0].NodeIndices(),
            (std::vector<fesa::EntityIndex>{4, 1, 6, 3, 7, 2, 5, 0}));
  auto model = fesa::AnalysisModel::Create(domain.Value());
  ASSERT_TRUE(model.HasValue());
  auto dofs = fesa::DofManager::Create(model.Value());
  ASSERT_TRUE(dofs.HasValue());
  EXPECT_EQ(dofs.Value().ConstrainedDofs(),
            (std::vector<std::size_t>{12, 13, 14}));
  EXPECT_FALSE(Map(Replace(deck, "20,1,0,0", "80,1,0,0")).HasValue());
  EXPECT_FALSE(
      Map(Replace(deck, "99,10,20,30,40,50,60,70,80",
                  "99,10,20,30,40,50,60,70,80\n99,10,20,30,40,50,60,70,80"))
          .HasValue());
  auto shuffled = SolidDefinition();
  std::swap(shuffled.nodes[0], shuffled.nodes[7]);
  shuffled.solid_elements[0].node_indices = {7, 1, 2, 3, 4, 5, 6, 0};
  auto direct = fesa::Domain::Create(std::move(shuffled));
  ASSERT_TRUE(direct.HasValue());
  EXPECT_EQ(direct.Value().Nodes()[0].source_id.source_label, 80);
  EXPECT_EQ(direct.Value().Nodes()[7].source_id.source_label, 10);
  EXPECT_EQ(direct.Value().SolidElements()[0].node_indices,
            (std::array<fesa::EntityIndex, 8>{7, 1, 2, 3, 4, 5, 6, 0}));
  auto nodes = SolidDefinition();
  nodes.nodes[1].source_id = nodes.nodes[0].source_id;
  EXPECT_FALSE(fesa::Domain::Create(std::move(nodes)).HasValue());
  auto elements = SolidDefinition();
  elements.solid_elements.push_back(elements.solid_elements[0]);
  EXPECT_FALSE(fesa::Domain::Create(std::move(elements)).HasValue());
}
TEST(Hex8Semantics, RejectsUnknownSolidAssignmentTargets) {
  for (const std::string assignment :
       {"elset=Missing,material=M", "elset=E,material=Missing"}) {
    SCOPED_TRACE(assignment);
    const auto result = Map(Replace(Deck(), "elset=E,material=M", assignment));
    ASSERT_FALSE(result.HasValue());
    ASSERT_FALSE(result.GetStatus().Diagnostics().empty());
    EXPECT_EQ(result.GetStatus().Diagnostics()[0].code,
              "unresolved-solid-section");
  }
}
TEST(Hex8Semantics, RejectsInvalidSolidDynamicSchedulesAndIncrementLimits) {
  const auto valid = Deck("*Dynamic,direct\n.1,1\n", "", "*Density\n2\n");
  ASSERT_TRUE(Map(valid).HasValue());
  for (const std::string schedule : {".3,1", "0,1", "-1,1", ".1,0", "nan,1"}) {
    SCOPED_TRACE(schedule);
    EXPECT_FALSE(Map(Replace(valid, ".1,1\n", schedule + "\n")).HasValue());
  }
  for (const std::string limit : {"0", "9", "4294967296"}) {
    SCOPED_TRACE(limit);
    EXPECT_FALSE(Map(Replace(valid, "inc=100", "inc=" + limit)).HasValue());
  }
}
}  // namespace