Skeleton for a fixed connectivity projection solver

src/scheme/FluxingAdvectionSolver.cpp

+#include <scheme/FluxingAdvectionSolver.hpp>
+
+#include <language/utils/EvaluateAtPoints.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/IMesh.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+template <size_t Dimension>
+class FluxingAdvectionSolver
+{
+ private:
+  using Rd = TinyVector<Dimension>;
+
+  using MeshType = Mesh<Connectivity<Dimension>>;
+
+  using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, double>;
+
+  const MeshType& m_old_mesh;
+  const std::shared_ptr<const MeshType> m_new_mesh;
+  const DiscreteFunctionP0<Dimension, const double> m_old_q;
+
+ public:
+  // CellValue<double>
+  // compute_PFnp1(const DiscreteFunctionP0<Dimension, const double> F, const double& dt, const double& dx)
+  // {
+  //   CellValue<double> PFnp1{m_mesh.connectivity()};
+
+  //   DiscreteFunctionP0<Dimension, double> deltaF  = compute_delta2Fn(F);
+  //   DiscreteFunctionP0<Dimension, double> deltaF0 = compute_delta2Fn(m_Fn);
+
+  //   for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+  //     PFnp1[cell_id] = m_Fn[cell_id][0] + m_Fn[cell_id][1] -
+  //                      (0.5 * dt / dx) * m_lambda * (deltaF[cell_id][0] - deltaF[cell_id][1]) -
+  //                      (0.5 * dt / dx) * m_lambda * (deltaF0[cell_id][0] - deltaF0[cell_id][1]);
+  //   }
+  //   return PFnp1;
+  // }
+
+  // DiscreteFunctionP0<Dimension, double>
+  // apply(const double& dt, const double& eps)
+  // {
+  //   const DiscreteFunctionP0<Dimension, const double>& F0 = m_Fn;
+  //   DiscreteFunctionP0<Dimension, double> Fnp1            = copy(F0);
+  //   DiscreteFunctionP0<Dimension, double> deltaFn         = compute_delta2Fn(F0);
+
+  //   for (size_t p = 0; p < 2; ++p) {
+  //     CellId first_cell_id = 0;
+  //     const double dx      = m_dx_table[first_cell_id];
+
+  //     DiscreteFunctionP0<Dimension, double> deltaFnp1 = compute_delta2Fn(Fnp1);
+
+  //     const CellValue<const double> PFnp1  = compute_PFnp1(Fnp1, dt, dx);
+  //     const CellValue<const double> APFnp1 = getA(PFnp1);
+  //     const CellArray<const double> MPFnp1 = compute_M(PFnp1, APFnp1);
+  //     const CellValue<const double> PFn    = compute_PFn(F0);
+  //     const CellValue<const double> APFn   = getA(PFn);
+  //     const CellArray<const double> MPFn   = compute_M(PFn, APFn);
+
+  //     for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+  //       Fnp1[cell_id][0] = 1. / (1 + 0.5 * dt / eps) *
+  //                          ((0.5 * dt / eps) * MPFnp1[cell_id][0] + F0[cell_id][0] -
+  //                           (0.5 * dt / dx) * m_lambda * (deltaFnp1[cell_id][0] + deltaFn[cell_id][0]) +
+  //                           (0.5 * dt / eps) * (MPFn[cell_id][0] - F0[cell_id][0]));
+  //       Fnp1[cell_id][1] = 1. / (1 + 0.5 * dt / eps) *
+  //                          ((0.5 * dt / eps) * MPFnp1[cell_id][1] + F0[cell_id][1] +
+  //                           (0.5 * dt / dx) * m_lambda * (deltaFnp1[cell_id][1] + deltaFn[cell_id][1]) +
+  //                           (0.5 * dt / eps) * (MPFn[cell_id][1] - F0[cell_id][1]));
+  //     }
+  //   }
+
+  //   return Fnp1;
+  // }
+
+  DiscreteFunctionP0<Dimension, double>
+  apply()
+  {
+    //    std::shared_ptr<const MeshType> mesh =
+    //  std::make_shared<MeshType>(m_new_mesh.shared_connectivity(), m_new_mesh.xr());
+    DiscreteFunctionP0<Dimension, double> new_q(m_new_mesh);
+    if (m_new_mesh->shared_connectivity() != m_old_mesh.shared_connectivity()) {
+      throw NormalError("Old and new meshes must share the same connectivity");
+    }
+    return new_q;
+  }
+
+  FluxingAdvectionSolver(const std::shared_ptr<const MeshType> old_mesh,
+                         const std::shared_ptr<const MeshType> new_mesh,
+                         const DiscreteFunctionP0<Dimension, const double>& old_q)
+    : m_old_mesh{*old_mesh}, m_new_mesh{new_mesh}, m_old_q{old_q}
+  {}
+
+  ~FluxingAdvectionSolver() = default;
+};
+
+std::shared_ptr<const IDiscreteFunction>
+FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
+                              const std::shared_ptr<const IDiscreteFunction>& old_q)
+{
+  const std::shared_ptr<const IMesh> old_mesh = getCommonMesh({old_q});
+  if (old_q->descriptor().type() != DiscreteFunctionType::P0) {
+    throw NormalError("Invalid discrete function type. Expecting P0.");
+  }
+
+  switch (old_mesh->dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+
+    const std::shared_ptr<const MeshType> old_mesh0 = std::dynamic_pointer_cast<const MeshType>(old_mesh);
+
+    const DiscreteFunctionP0<Dimension, double>& old_q0 =
+      dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*old_q);
+
+    const std::shared_ptr<const MeshType> new_mesh0 = std::dynamic_pointer_cast<const MeshType>(new_mesh);
+
+    FluxingAdvectionSolver<Dimension> solver(old_mesh0, new_mesh0, old_q0);
+
+    return std::make_shared<const DiscreteFunctionP0<Dimension, double>>(solver.apply());
+  }
+  case 2: {
+    throw NotImplementedError("Fluxing advection solver not implemented in dimension 2");
+  }
+  case 3: {
+    throw NotImplementedError("Fluxing advection solver not implemented in dimension 3");
+  }
+  default: {
+    throw UnexpectedError("Invalid mesh dimension");
+  }
+  }
+}

src/scheme/FluxingAdvectionSolver.hpp

+#ifndef FLUXING_ADVECION_SOLVER_HPP
+#define FLUXING_ADVECION_SOLVER_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+
+std::shared_ptr<const IDiscreteFunction> FluxingAdvectionSolverHandler(
+  const std::shared_ptr<const IMesh> new_mesh,
+  const std::shared_ptr<const IDiscreteFunction>& old_q);
+
+#endif   // FLUXING_ADVECION_SOLVER_HPP