#include <meson_mfem_template/config.hpp>
#include <mfem.hpp>

#include <cmath>
#include <iostream>

int main(int argc, char **argv)
{
   const char *device_name = argc > 1 ? argv[1] :
#if MESON_MFEM_HAS_CEED && !defined(__EMSCRIPTEN__)
                             "ceed-cpu";
#else
                             "cpu";
#endif
   mfem::Device device(device_name);

   mfem::Mesh mesh = mfem::Mesh::MakeCartesian2D(
      6, 6, mfem::Element::QUADRILATERAL, true, 1.0, 1.0);
   const int dimension = mesh.Dimension();
   mfem::H1_FECollection elements(2, dimension);
   mfem::FiniteElementSpace space(&mesh, &elements);

   mfem::Array<int> essential_boundary(mesh.bdr_attributes.Max());
   essential_boundary = 1;
   mfem::Array<int> essential_dofs;
   space.GetEssentialTrueDofs(essential_boundary, essential_dofs);

   mfem::ConstantCoefficient one(1.0);
   mfem::LinearForm rhs(&space);
   rhs.AddDomainIntegrator(new mfem::DomainLFIntegrator(one));
   rhs.Assemble();

   mfem::GridFunction solution(&space);
   solution = 0.0;

   mfem::BilinearForm diffusion(&space);
   diffusion.AddDomainIntegrator(new mfem::DiffusionIntegrator(one));
   if (mfem::Device::Allows(mfem::Backend::CEED_MASK))
   {
      diffusion.SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
   }
   diffusion.Assemble();

   mfem::OperatorPtr matrix;
   mfem::Vector linear_rhs;
   mfem::Vector linear_solution;
   diffusion.FormLinearSystem(
      essential_dofs, solution, rhs, matrix, linear_solution, linear_rhs);

   mfem::CGSolver solver;
   solver.SetOperator(*matrix);
   solver.SetRelTol(1e-12);
   solver.SetAbsTol(0.0);
   solver.SetMaxIter(200);
   solver.SetPrintLevel(0);
   solver.Mult(linear_rhs, linear_solution);
   diffusion.RecoverFEMSolution(linear_solution, rhs, solution);

#if MESON_MFEM_HAS_FMS
   mfem::FMSDataCollection fms_collection("meson-mfem-smoke", &mesh);
   fms_collection.SetProtocol("ascii");
#endif

#if MESON_MFEM_HAS_ALGOIM
   mfem::FunctionCoefficient level_set(
      [](const mfem::Vector &point) { return point[0] + point[1] - 0.2; });
   mfem::AlgoimIntegrationRules algoim_rules(2, level_set, 2);
   mfem::IntegrationRule cut_rule;
   algoim_rules.GetVolumeIntegrationRule(
      *mesh.GetElementTransformation(0), cut_rule);
   if (cut_rule.GetNPoints() == 0)
   {
      return 3;
   }
#endif

   const double norm = solution.Norml2();
   std::cout << "MFEM serial Poisson: dofs=" << space.GetTrueVSize()
             << " l2=" << norm << " device=" << device_name << '\n';
   return solver.GetConverged() && std::isfinite(norm) && norm > 0.0 ? 0 : 2;
}