#include <PugsUtils.hpp>

#include <rang.hpp>

#include <Connectivity.hpp>

#include <AcousticSolver.hpp>
#include <BoundaryCondition.hpp>
#include <Mesh.hpp>

#include <VTKWriter.hpp>

#include <Timer.hpp>

#include <TinyMatrix.hpp>
#include <TinyVector.hpp>

#include <BoundaryConditionDescriptor.hpp>

#include <MeshNodeBoundary.hpp>

#include <GmshReader.hpp>
#include <PugsParser.hpp>

#include <SynchronizerManager.hpp>

#include <iostream>

#include <limits>
#include <map>
#include <regex>

int
main(int argc, char* argv[])
{
  std::string filename = initialize(argc, argv);

  std::regex gmsh_regex("(.*).msh");
  if (not std::regex_match(filename, gmsh_regex)) {
    SynchronizerManager::create();
    parser(filename);
    SynchronizerManager::destroy();
    return 0;
  }

  std::map<std::string, double> method_cost_map;

  SynchronizerManager::create();

  if (filename != "") {
    std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
    Timer gmsh_timer;
    gmsh_timer.reset();
    GmshReader gmsh_reader(filename);
    method_cost_map["Mesh building"] = gmsh_timer.seconds();

    std::shared_ptr<IMesh> p_mesh = gmsh_reader.mesh();

    switch (p_mesh->dimension()) {
    case 1: {
      std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
      std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
      for (const auto& sym_boundary_name : sym_boundary_name_list) {
        std::shared_ptr<BoundaryDescriptor> boudary_descriptor =
          std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
        SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
          new SymmetryBoundaryConditionDescriptor(boudary_descriptor);

        bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
      }

      using ConnectivityType = Connectivity1D;
      using MeshType         = Mesh<ConnectivityType>;
      using MeshDataType     = MeshData<MeshType>;
      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;

      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());

      Timer timer;
      timer.reset();
      MeshDataType mesh_data(mesh);

      std::vector<BoundaryConditionHandler> bc_list;
      {
        for (const auto& bc_descriptor : bc_descriptor_list) {
          switch (bc_descriptor->type()) {
          case BoundaryConditionDescriptor::Type::symmetry: {
            const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
              dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
            for (size_t i_ref_node_list = 0;
                 i_ref_node_list < mesh.connectivity().numberOfRefItemList<ItemType::node>(); ++i_ref_node_list) {
              const RefNodeList& ref_node_list = mesh.connectivity().refItemList<ItemType::node>(i_ref_node_list);
              const RefId& ref                 = ref_node_list.refId();
              if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
                  new SymmetryBoundaryCondition<MeshType::Dimension>(
                    MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_node_list));
                std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
                bc_list.push_back(BoundaryConditionHandler(bc));
              }
            }
            break;
          }
          default: {
            std::cerr << "Unknown BCDescription\n";
            std::exit(1);
          }
          }
        }
      }

      UnknownsType unknowns(mesh_data);

      unknowns.initializeSod();

      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);

      using Rd = TinyVector<MeshType::Dimension>;

      const CellValue<const double>& Vj = mesh_data.Vj();

      const double tmax = 0.2;
      double t          = 0;

      int itermax   = std::numeric_limits<int>::max();
      int iteration = 0;

      CellValue<double>& rhoj   = unknowns.rhoj();
      CellValue<double>& ej     = unknowns.ej();
      CellValue<double>& pj     = unknowns.pj();
      CellValue<double>& gammaj = unknowns.gammaj();
      CellValue<double>& cj     = unknowns.cj();

      BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);

      VTKWriter vtk_writer("mesh", 0.01);

      while ((t < tmax) and (iteration < itermax)) {
        vtk_writer.write(mesh,
                         {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                          NamedItemValue{"coords", mesh.xr()},
                          NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                          NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
                         t);
        double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
        if (t + dt > tmax) {
          dt = tmax - t;
        }
        acoustic_solver.computeNextStep(t, dt, unknowns);

        block_eos.updatePandCFromRhoE();

        t += dt;
        ++iteration;
      }
      vtk_writer.write(mesh,
                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                        NamedItemValue{"coords", mesh.xr()},
                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
                       t, true);   // forces last output

      std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
                << t << rang::fg::reset << " (" << iteration << " iterations)\n";

      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();

      {   // gnuplot output for density
        const CellValue<const Rd>& xj       = mesh_data.xj();
        const CellValue<const double>& rhoj = unknowns.rhoj();
        std::ofstream fout("rho");
        for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
          fout << xj[j][0] << ' ' << rhoj[j] << '\n';
        }
      }

      break;
    }
    case 2: {
      // test case boundary condition description
      std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX"};
      std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
      for (const auto& sym_boundary_name : sym_boundary_name_list) {
        std::shared_ptr<BoundaryDescriptor> boudary_descriptor =
          std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
        SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
          new SymmetryBoundaryConditionDescriptor(boudary_descriptor);

        bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
      }

      using ConnectivityType = Connectivity2D;
      using MeshType         = Mesh<ConnectivityType>;
      using MeshDataType     = MeshData<MeshType>;
      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;

      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());

      Timer timer;
      timer.reset();
      MeshDataType mesh_data(mesh);

      std::vector<BoundaryConditionHandler> bc_list;
      {
        for (const auto& bc_descriptor : bc_descriptor_list) {
          switch (bc_descriptor->type()) {
          case BoundaryConditionDescriptor::Type::symmetry: {
            const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
              dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
            for (size_t i_ref_face_list = 0;
                 i_ref_face_list < mesh.connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
              const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
              const RefId& ref                 = ref_face_list.refId();
              if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
                  new SymmetryBoundaryCondition<MeshType::Dimension>(
                    MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
                std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
                bc_list.push_back(BoundaryConditionHandler(bc));
              }
            }
            break;
          }
          default: {
            std::cerr << "Unknown BCDescription\n";
            std::exit(1);
          }
          }
        }
      }

      UnknownsType unknowns(mesh_data);

      unknowns.initializeSod();

      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);

      const CellValue<const double>& Vj = mesh_data.Vj();

      const double tmax = 0.2;
      double t          = 0;

      int itermax   = std::numeric_limits<int>::max();
      int iteration = 0;

      CellValue<double>& rhoj   = unknowns.rhoj();
      CellValue<double>& ej     = unknowns.ej();
      CellValue<double>& pj     = unknowns.pj();
      CellValue<double>& gammaj = unknowns.gammaj();
      CellValue<double>& cj     = unknowns.cj();

      BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);

      VTKWriter vtk_writer("mesh", 0.01);

      while ((t < tmax) and (iteration < itermax)) {
        vtk_writer.write(mesh,
                         {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                          NamedItemValue{"coords", mesh.xr()},
                          NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                          NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
                         t);
        double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
        if (t + dt > tmax) {
          dt = tmax - t;
        }
        acoustic_solver.computeNextStep(t, dt, unknowns);

        block_eos.updatePandCFromRhoE();

        t += dt;
        ++iteration;
      }
      vtk_writer.write(mesh,
                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                        NamedItemValue{"coords", mesh.xr()},
                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
                       t, true);   // forces last output

      std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
                << t << rang::fg::reset << " (" << iteration << " iterations)\n";

      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
      break;
    }
    case 3: {
      std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
      std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
      for (const auto& sym_boundary_name : sym_boundary_name_list) {
        std::shared_ptr<BoundaryDescriptor> boudary_descriptor =
          std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
        SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
          new SymmetryBoundaryConditionDescriptor(boudary_descriptor);

        bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
      }

      using ConnectivityType = Connectivity3D;
      using MeshType         = Mesh<ConnectivityType>;
      using MeshDataType     = MeshData<MeshType>;
      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;

      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());

      Timer timer;
      timer.reset();
      MeshDataType mesh_data(mesh);

      std::vector<BoundaryConditionHandler> bc_list;
      {
        for (const auto& bc_descriptor : bc_descriptor_list) {
          switch (bc_descriptor->type()) {
          case BoundaryConditionDescriptor::Type::symmetry: {
            const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
              dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
            for (size_t i_ref_face_list = 0;
                 i_ref_face_list < mesh.connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
              const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
              const RefId& ref                 = ref_face_list.refId();
              if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
                  new SymmetryBoundaryCondition<MeshType::Dimension>(
                    MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
                std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
                bc_list.push_back(BoundaryConditionHandler(bc));
              }
            }
            break;
          }
          default: {
            std::cerr << "Unknown BCDescription\n";
            std::exit(1);
          }
          }
        }
      }

      UnknownsType unknowns(mesh_data);

      unknowns.initializeSod();

      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);

      const CellValue<const double>& Vj = mesh_data.Vj();

      const double tmax = 0.2;
      double t          = 0;

      int itermax   = std::numeric_limits<int>::max();
      int iteration = 0;

      CellValue<double>& rhoj   = unknowns.rhoj();
      CellValue<double>& ej     = unknowns.ej();
      CellValue<double>& pj     = unknowns.pj();
      CellValue<double>& gammaj = unknowns.gammaj();
      CellValue<double>& cj     = unknowns.cj();

      BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);

      VTKWriter vtk_writer("mesh", 0.01);

      while ((t < tmax) and (iteration < itermax)) {
        vtk_writer.write(mesh,
                         {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                          NamedItemValue{"coords", mesh.xr()},
                          NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                          NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
                         t);
        double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
        if (t + dt > tmax) {
          dt = tmax - t;
        }
        acoustic_solver.computeNextStep(t, dt, unknowns);
        block_eos.updatePandCFromRhoE();

        t += dt;
        ++iteration;
      }
      vtk_writer.write(mesh,
                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                        NamedItemValue{"coords", mesh.xr()},
                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
                       t, true);   // forces last output

      std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
                << t << rang::fg::reset << " (" << iteration << " iterations)\n";

      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
      break;
    }
    }

    std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':'
              << __LINE__ << ")\n";

  } else {
    std::cerr << "Connectivity1D defined by number of nodes no more implemented\n";
    std::exit(0);
  }

  SynchronizerManager::destroy();

  finalize();

  std::string::size_type size = 0;
  for (const auto& method_cost : method_cost_map) {
    size = std::max(size, method_cost.first.size());
  }

  for (const auto& method_cost : method_cost_map) {
    std::cout << "* [" << rang::fgB::cyan << std::setw(size) << std::left << method_cost.first << rang::fg::reset
              << "] Execution time: " << rang::style::bold << method_cost.second << rang::style::reset << '\n';
  }

  return 0;
}