From 26a054156f83a1d4df0bc460ca11cef8d3ae9860 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 20 Mar 2023 22:56:00 +0100
Subject: [PATCH] Improve fluxing based remapping

- remap tuples of Vh
- handle remapping for P0 vectors
- add remapping for 1d
---
 src/language/modules/SchemeModule.cpp |  27 +-
 src/scheme/FluxingAdvectionSolver.cpp | 616 +++++++++++++++-----------
 src/scheme/FluxingAdvectionSolver.hpp |   6 +-
 3 files changed, 356 insertions(+), 293 deletions(-)

diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 8363afeb0..cf183409d 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -615,23 +615,16 @@ SchemeModule::SchemeModule()
                                                  }
 
                                                  ));
-  this
-    ->_addBuiltinFunction("fluxing_advection",
-                          std::function(
-
-                            // [](const std::shared_ptr<const IMesh> new_mesh,
-                            //    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& remapped_variables)
-                            //   -> std::vector<std::shared_ptr<const DiscreteFunctionVariant>> {
-                            //   return FluxingAdvectionSolverHandler(new_mesh, remapped_variables);
-                            // }
-
-                            [](const std::shared_ptr<const IMesh> new_mesh,
-                               const std::shared_ptr<const DiscreteFunctionVariant>& remapped_variables)
-                              -> std::shared_ptr<const DiscreteFunctionVariant> {
-                              return FluxingAdvectionSolverHandler(new_mesh, remapped_variables);
-                            }
-
-                            ));
+  this->_addBuiltinFunction("fluxing_advection",
+                            std::function(
+
+                              [](const std::shared_ptr<const IMesh> new_mesh,
+                                 const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& remapped_variables)
+                                -> std::vector<std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return advectByFluxing(new_mesh, remapped_variables);
+                              }
+
+                              ));
 
   MathFunctionRegisterForVh{*this};
 }
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index c4884b62b..d2f7cd45c 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -11,6 +11,7 @@
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
+
 template <size_t Dimension>
 class FluxingAdvectionSolver
 {
@@ -23,359 +24,432 @@ class FluxingAdvectionSolver
   const std::shared_ptr<const MeshType> m_old_mesh;
   const std::shared_ptr<const MeshType> m_new_mesh;
 
+  using RemapVariant = std::variant<CellValue<double>,
+                                    CellValue<TinyVector<1>>,
+                                    CellValue<TinyVector<2>>,
+                                    CellValue<TinyVector<3>>,
+                                    CellValue<TinyMatrix<1>>,
+                                    CellValue<TinyMatrix<2>>,
+                                    CellValue<TinyMatrix<3>>,
+
+                                    CellArray<double>>;
+
+  std::vector<RemapVariant> m_remapped_list;
+
+  FaceValue<const CellId> m_donnor_cell;
+  FaceValue<const double> m_cycle_fluxing_volume;
+  size_t m_number_of_cycles;
+
+  FaceValue<double> _computeAlgebraicFluxingVolume();
+  void _computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes);
+  FaceValue<double> _computeFluxingVolume(FaceValue<double> algebraic_fluxing_volumes);
+  void _computeCycleNumber(FaceValue<double> fluxing_volumes);
+  void _computeGeometricalData();
+
+  template <typename DataType>
+  void
+  _storeValues(const DiscreteFunctionP0<Dimension, const DataType>& old_q)
+  {
+    m_remapped_list.emplace_back(copy(old_q.cellValues()));
+  }
+
+  template <typename DataType>
+  void
+  _storeValues(const DiscreteFunctionP0Vector<Dimension, const DataType>& old_q)
+  {
+    m_remapped_list.emplace_back(copy(old_q.cellArrays()));
+  }
+
+  template <typename CellDataType>
+  void _remapOne(const CellValue<const double>& step_Vj, CellDataType& old_q);
+
+  void _remapAllQuantities();
+
  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;
-  // }
-
-  FaceValue<double> computeFluxVolume() const;
-
-  FluxingAdvectionSolver(const std::shared_ptr<const MeshType> old_mesh, const std::shared_ptr<const MeshType> new_mesh)
-    : m_old_mesh{old_mesh}, m_new_mesh{new_mesh}
-  {}
+  std::vector<std::shared_ptr<const DiscreteFunctionVariant>>   //
+  remap(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& quantities);
+
+  FluxingAdvectionSolver(const std::shared_ptr<const IMesh> i_old_mesh, const std::shared_ptr<const IMesh> i_new_mesh)
+    : m_old_mesh{std::dynamic_pointer_cast<const MeshType>(i_old_mesh)},
+      m_new_mesh{std::dynamic_pointer_cast<const MeshType>(i_new_mesh)}
+  {
+    if ((m_old_mesh.use_count() == 0) or (m_new_mesh.use_count() == 0)) {
+      throw NormalError("old and new meshes must be of same type");
+    }
+
+    if (m_new_mesh->shared_connectivity() != m_old_mesh->shared_connectivity()) {
+      throw NormalError("old and new meshes must share the same connectivity");
+    }
+
+    this->_computeGeometricalData();
+  }
 
   ~FluxingAdvectionSolver() = default;
 };
 
+template <size_t Dimension>
+void
+FluxingAdvectionSolver<Dimension>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
+{
+  m_donnor_cell = [&] {
+    const FaceValuePerCell<const bool> cell_face_is_reversed = m_new_mesh->connectivity().cellFaceIsReversed();
+    const auto face_to_cell_matrix                           = m_new_mesh->connectivity().faceToCellMatrix();
+
+    const auto face_local_number_in_their_cells = m_new_mesh->connectivity().faceLocalNumbersInTheirCells();
+
+    FaceValue<CellId> donnor_cell{m_old_mesh->connectivity()};
+    parallel_for(
+      m_new_mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+        const auto& face_to_cell = face_to_cell_matrix[face_id];
+        if (face_to_cell.size() == 1) {
+          donnor_cell[face_id] = face_to_cell[0];
+        } else {
+          const CellId cell_id        = face_to_cell[0];
+          const size_t i_face_in_cell = face_local_number_in_their_cells[face_id][0];
+          if (cell_face_is_reversed[cell_id][i_face_in_cell] xor (algebraic_fluxing_volumes[face_id] <= 0)) {
+            donnor_cell[face_id] = cell_id;
+          } else {
+            donnor_cell[face_id] = face_to_cell[1];
+          }
+        }
+      });
+
+    return donnor_cell;
+  }();
+}
+
+template <>
+void
+FluxingAdvectionSolver<1>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
+{
+  m_donnor_cell = [&] {
+    const auto face_to_cell_matrix = m_new_mesh->connectivity().faceToCellMatrix();
+    const auto cell_to_face_matrix = m_new_mesh->connectivity().cellToFaceMatrix();
+
+    FaceValue<CellId> donnor_cell{m_old_mesh->connectivity()};
+    parallel_for(
+      m_new_mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+        const auto& face_to_cell = face_to_cell_matrix[face_id];
+        if (face_to_cell.size() == 1) {
+          donnor_cell[face_id] = face_to_cell[0];
+        } else {
+          const CellId cell_id = face_to_cell[0];
+          if ((algebraic_fluxing_volumes[face_id] <= 0) xor (cell_to_face_matrix[cell_id][0] == face_id)) {
+            donnor_cell[face_id] = cell_id;
+          } else {
+            donnor_cell[face_id] = face_to_cell[1];
+          }
+        }
+      });
+
+    return donnor_cell;
+  }();
+}
+
 template <>
 FaceValue<double>
-FluxingAdvectionSolver<1>::computeFluxVolume() const
+FluxingAdvectionSolver<1>::_computeAlgebraicFluxingVolume()
 {
-  throw NotImplementedError("Viens");
+  Array<double> algebraic_fluxing_volumes{m_new_mesh->numberOfNodes()};
+  NodeValue<double> fluxing_volume(m_new_mesh->connectivity(), algebraic_fluxing_volumes);
+  auto old_xr = m_old_mesh->xr();
+  auto new_xr = m_new_mesh->xr();
+
+  parallel_for(
+    m_new_mesh->numberOfNodes(),
+    PUGS_LAMBDA(NodeId node_id) { fluxing_volume[node_id] = new_xr[node_id][0] - old_xr[node_id][0]; });
+
+  return FaceValue<double>{m_new_mesh->connectivity(), algebraic_fluxing_volumes};
 }
 
 template <>
 FaceValue<double>
-FluxingAdvectionSolver<2>::computeFluxVolume() const
+FluxingAdvectionSolver<2>::_computeAlgebraicFluxingVolume()
 {
-  if (m_new_mesh->shared_connectivity() != m_old_mesh->shared_connectivity()) {
-    throw NormalError("Old and new meshes must share the same connectivity");
-  }
-  // std::cout << " CARRE "
-  // << "\n";
-  // MeshDataType& old_mesh_data    = MeshDataManager::instance().getMeshData(*m_old_mesh);
-  // MeshDataType& new_mesh_data    = MeshDataManager::instance().getMeshData(*m_new_mesh);
   const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix();
-  FaceValue<double> flux_volume(m_new_mesh->connectivity());
+  FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity());
+  auto old_xr = m_old_mesh->xr();
+  auto new_xr = m_new_mesh->xr();
   parallel_for(
     m_new_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-      const auto face_to_node = face_to_node_matrix[face_id];
-      const Rd x0             = m_old_mesh->xr()[face_to_node[0]];
-      const Rd x1             = m_old_mesh->xr()[face_to_node[1]];
-      const Rd x2             = m_new_mesh->xr()[face_to_node[1]];
-      const Rd x3             = m_new_mesh->xr()[face_to_node[0]];
-      TinyMatrix<2> M(x2[0] - x0[0], x3[0] - x1[0], x2[1] - x0[1], x3[1] - x1[1]);
-      flux_volume[face_id] = 0.5 * det(M);
-      // std::cout << " x1 " << x1 << " x0 " << x0 << " x3 " << x3 << " x2 " << x2 << " flux volume "
-      //           << flux_volume[face_id] << "\n";
+      const auto& face_to_node = face_to_node_matrix[face_id];
+
+      const Rd& x0 = old_xr[face_to_node[0]];
+      const Rd& x1 = old_xr[face_to_node[1]];
+      const Rd& x2 = new_xr[face_to_node[1]];
+      const Rd& x3 = new_xr[face_to_node[0]];
+
+      TinyMatrix<2> M(x3[0] - x1[0], x2[0] - x0[0],   //
+                      x3[1] - x1[1], x2[1] - x0[1]);
+
+      algebraic_fluxing_volume[face_id] = 0.5 * det(M);
     });
-  return flux_volume;
+
+  return algebraic_fluxing_volume;
 }
 
 template <>
 FaceValue<double>
-FluxingAdvectionSolver<3>::computeFluxVolume() const
+FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
 {
   throw NotImplementedError("ViensViensViens");
 }
 
-template <typename MeshType, typename DataType>
-auto
-calculateRemapCycles(const std::shared_ptr<const MeshType>& old_mesh,
-                     [[maybe_unused]] const FaceValue<DataType>& fluxing_volumes)
+template <size_t Dimension>
+FaceValue<double>
+FluxingAdvectionSolver<Dimension>::_computeFluxingVolume(FaceValue<double> algebraic_fluxing_volumes)
+{
+  Assert(m_donnor_cell.isBuilt());
+  // Now that donnor cells are clearly defined, we consider the
+  // non-algebraic volumes of fluxing
+  parallel_for(
+    algebraic_fluxing_volumes.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) {
+      algebraic_fluxing_volumes[face_id] = std::abs(algebraic_fluxing_volumes[face_id]);
+    });
+
+  return algebraic_fluxing_volumes;
+}
+
+template <size_t Dimension>
+void
+FluxingAdvectionSolver<Dimension>::_computeCycleNumber(FaceValue<double> fluxing_volumes)
 {
-  constexpr size_t Dimension                               = MeshType::Dimension;
-  const FaceValuePerCell<const bool> cell_face_is_reversed = old_mesh->connectivity().cellFaceIsReversed();
-  const auto cell_to_face_matrix                           = old_mesh->connectivity().cellToFaceMatrix();
-  const CellValue<double> total_negative_flux(old_mesh->connectivity());
+  const auto cell_to_face_matrix = m_old_mesh->connectivity().cellToFaceMatrix();
+
+  const CellValue<double> total_negative_flux(m_old_mesh->connectivity());
   total_negative_flux.fill(0);
+
   parallel_for(
-    old_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+    m_old_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
       const auto& cell_to_face = cell_to_face_matrix[cell_id];
       for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
         FaceId face_id = cell_to_face[i_face];
-        double flux    = fluxing_volumes[face_id];
-        if (!cell_face_is_reversed(cell_id, i_face)) {
-          flux = -flux;
-        }
-        if (flux < 0) {
-          total_negative_flux[cell_id] += flux;
+        if (cell_id == m_donnor_cell[face_id]) {
+          total_negative_flux[cell_id] += fluxing_volumes[face_id];
         }
       }
-      // std::cout << " cell_id " << cell_id << " total_negative_flux " << total_negative_flux[cell_id] << "\n";
     });
-  MeshData<Dimension>& mesh_data   = MeshDataManager::instance().getMeshData(*old_mesh);
+
+  MeshData<Dimension>& mesh_data   = MeshDataManager::instance().getMeshData(*m_old_mesh);
   const CellValue<const double> Vj = mesh_data.Vj();
-  const CellValue<size_t> ratio(old_mesh->connectivity());
+  const CellValue<size_t> ratio(m_old_mesh->connectivity());
+
   parallel_for(
-    old_mesh->numberOfCells(),
-    PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = std::ceil(std::abs(total_negative_flux[cell_id]) / Vj[cell_id]); });
+    m_old_mesh->numberOfCells(),
+    PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = std::ceil(total_negative_flux[cell_id] / Vj[cell_id]); });
 
-  return max(ratio);
-}
+  size_t number_of_cycles = max(ratio);
+
+  if (number_of_cycles > 1) {
+    const double cycle_ratio = 1. / number_of_cycles;
 
-template <typename MeshType, typename DataType>
-auto
-remapUsingFluxing(const std::shared_ptr<const MeshType>& old_mesh,
-                  const std::shared_ptr<const MeshType>& new_mesh,
-                  const FaceValue<double>& fluxing_volumes,
-                  const size_t num,
-                  const DiscreteFunctionP0<MeshType::Dimension, const DataType>& old_q)
-{
-  constexpr size_t Dimension = MeshType::Dimension;
-  //  const Connectivity<Dimension>& connectivity = new_mesh->connectivity();
-  const FaceValuePerCell<const bool> cell_face_is_reversed = new_mesh->connectivity().cellFaceIsReversed();
-  DiscreteFunctionP0<Dimension, DataType> new_q(new_mesh, copy(old_q.cellValues()));
-  DiscreteFunctionP0<Dimension, DataType> previous_q(new_mesh, copy(old_q.cellValues()));
-  const auto cell_to_face_matrix      = new_mesh->connectivity().cellToFaceMatrix();
-  const auto face_to_cell_matrix      = new_mesh->connectivity().faceToCellMatrix();
-  MeshData<Dimension>& old_mesh_data  = MeshDataManager::instance().getMeshData(*old_mesh);
-  const CellValue<const double> oldVj = old_mesh_data.Vj();
-  const CellValue<double> Vjstep(new_mesh->connectivity());
-  parallel_for(
-    new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-      Vjstep[cell_id] = oldVj[cell_id];
-      new_q[cell_id] *= oldVj[cell_id];
-    });
-  for (size_t jstep = 0; jstep < num; ++jstep) {
-    // std::cout << " step " << jstep << "\n";
-    parallel_for(
-      new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-        const auto& cell_to_face = cell_to_face_matrix[cell_id];
-        for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
-          FaceId face_id = cell_to_face[i_face];
-          double flux    = fluxing_volumes[face_id];
-          if (!cell_face_is_reversed(cell_id, i_face)) {
-            flux = -flux;
-          }
-          const auto& face_to_cell = face_to_cell_matrix[face_id];
-          if (face_to_cell.size() == 1) {
-            continue;
-          }
-          CellId other_cell_id = face_to_cell[0];
-          if (other_cell_id == cell_id) {
-            other_cell_id = face_to_cell[1];
-          }
-          DataType fluxed_q = previous_q[cell_id];
-          if (flux > 0) {
-            fluxed_q = previous_q[other_cell_id];
-          }
-          Vjstep[cell_id] += flux / num;
-          fluxed_q *= flux / num;
-          new_q[cell_id] += fluxed_q;
-        }
-        // std::cout << " old q " << old_q[cell_id] << " new q " << new_q[cell_id] << "\n";
-      });
     parallel_for(
-      new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-        previous_q[cell_id] = 1 / Vjstep[cell_id] * new_q[cell_id];
-        //     std::cout << " old q " << old_q[cell_id] << " new q " << previous_q[cell_id] << "\n";
-        //     std::cout << " old vj " << oldVj[cell_id] << " new Vj " << Vjstep[cell_id] << "\n";
-      });
+      fluxing_volumes.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) { fluxing_volumes[face_id] *= cycle_ratio; });
   }
 
-  MeshData<Dimension>& new_mesh_data  = MeshDataManager::instance().getMeshData(*new_mesh);
-  const CellValue<const double> newVj = new_mesh_data.Vj();
-  parallel_for(
-    new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { new_q[cell_id] = 1 / newVj[cell_id] * new_q[cell_id]; });
-  // for (CellId cell_id = 0; cell_id < new_mesh->numberOfCells(); ++cell_id) {
-  //   if (abs(newVj[cell_id] - Vjstep[cell_id]) > 1e-15) {
-  //     std::cout << " cell " << cell_id << " newVj " << newVj[cell_id] << " Vjstep " << Vjstep[cell_id] << " diff "
-  //               << abs(newVj[cell_id] - Vjstep[cell_id]) << "\n";
-  //   }
-  // }
-  return new_q;
+  m_number_of_cycles = number_of_cycles;
+  std::cout << " number_of_cycle " << number_of_cycles << "\n";
+
+  m_cycle_fluxing_volume = fluxing_volumes;
 }
 
-template <typename MeshType, typename DataType>
-auto
-remapUsingFluxing([[maybe_unused]] const std::shared_ptr<const MeshType>& old_mesh,
-                  const std::shared_ptr<const MeshType>& new_mesh,
-                  [[maybe_unused]] const FaceValue<double>& fluxing_volumes,
-                  [[maybe_unused]] const size_t num,
-                  const DiscreteFunctionP0Vector<MeshType::Dimension, const DataType>& old_q)
+template <size_t Dimension>
+void
+FluxingAdvectionSolver<Dimension>::_computeGeometricalData()
 {
-  constexpr size_t Dimension = MeshType::Dimension;
-  //  const Connectivity<Dimension>& connectivity = new_mesh->connectivity();
-
-  DiscreteFunctionP0Vector<Dimension, DataType> new_q(new_mesh, copy(old_q.cellArrays()));
-
-  throw NotImplementedError("DiscreteFunctionP0Vector");
-
-  return new_q;
+  auto fluxing_volumes = this->_computeAlgebraicFluxingVolume();
+  this->_computeDonorCells(fluxing_volumes);
+  fluxing_volumes = this->_computeFluxingVolume(fluxing_volumes);
+  this->_computeCycleNumber(fluxing_volumes);
 }
 
-std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
-FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
-                              const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& remapped_variables)
+template <size_t Dimension>
+template <typename CellDataType>
+void
+FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step_Vj, CellDataType& old_q)
 {
-  const std::shared_ptr<const IMesh> old_mesh = getCommonMesh(remapped_variables);
+  static_assert(is_item_value_v<CellDataType> or is_item_array_v<CellDataType>, "invalid data type");
 
-  if (not checkDiscretizationType({remapped_variables}, DiscreteFunctionType::P0)) {
-    throw NormalError("acoustic solver expects P0 functions");
-  }
+  const auto cell_to_face_matrix = m_new_mesh->connectivity().cellToFaceMatrix();
+  const auto face_to_cell_matrix = m_new_mesh->connectivity().faceToCellMatrix();
 
-  switch (old_mesh->dimension()) {
-  case 1: {
-    constexpr size_t Dimension = 1;
-    using MeshType             = Mesh<Connectivity<Dimension>>;
+  auto new_q = copy(old_q);
 
-    const std::shared_ptr<const MeshType> old_mesh0 = std::dynamic_pointer_cast<const MeshType>(old_mesh);
+  if constexpr (is_item_value_v<CellDataType>) {
+    parallel_for(
+      m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { new_q[cell_id] *= step_Vj[cell_id]; });
+  } else if constexpr (is_item_array_v<CellDataType>) {
+    parallel_for(
+      m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        auto new_array      = new_q[cell_id];
+        const double volume = step_Vj[cell_id];
 
-    const std::shared_ptr<const MeshType> new_mesh0 = std::dynamic_pointer_cast<const MeshType>(new_mesh);
+        for (size_t i = 0; i < new_array.size(); ++i) {
+          new_array[i] *= volume;
+        }
+      });
+  }
 
-    FluxingAdvectionSolver<Dimension> solver(old_mesh0, new_mesh0);
+  parallel_for(
+    m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+      const auto& cell_to_face = cell_to_face_matrix[cell_id];
+      for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+        const FaceId face_id        = cell_to_face[i_face];
+        const double fluxing_volume = m_cycle_fluxing_volume[face_id];
 
-    FaceValue<double> fluxing_volumes(new_mesh0->connectivity());
-    fluxing_volumes.fill(0);
+        const auto& face_to_cell = face_to_cell_matrix[face_id];
 
-    std::vector<std::shared_ptr<const DiscreteFunctionVariant>> new_variables;
+        if (face_to_cell.size() == 1) {
+          continue;
+        }
 
-    // for (auto&& variable_v : remapped_variables) {
-    //   std::visit(
-    //     [&](auto&& variable) {
-    //       using DiscreteFunctionT = std::decay_t<decltype(variable)>;
-    //       if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
-    //         remapUsingFluxing(new_mesh0, fluxing_volumes, variable);
-    //       }
-    //     },
-    //     variable_v->discreteFunction());
-    // }
+        CellId donnor_id = m_donnor_cell[face_id];
 
-    return remapped_variables;   // std::make_shared<std::vector<std::shared_ptr<const
-                                 // DiscreteFunctionVariant>>>(new_variables);
-  }
-  case 2: {
-    constexpr size_t Dimension = 2;
+        if constexpr (is_item_value_v<CellDataType>) {
+          auto fluxed_q = old_q[donnor_id];
+          fluxed_q *= ((donnor_id == cell_id) ? -1 : 1) * fluxing_volume;
 
-    using MeshType = Mesh<Connectivity<Dimension>>;
+          new_q[cell_id] += fluxed_q;
+        } else if constexpr (is_item_array_v<CellDataType>) {
+          const double sign   = ((donnor_id == cell_id) ? -1 : 1);
+          auto old_cell_array = old_q[donnor_id];
+          auto new_cell_array = new_q[cell_id];
+          for (size_t i = 0; i < new_cell_array.size(); ++i) {
+            new_cell_array[i] += (sign * fluxing_volume) * old_cell_array[i];
+          }
+        }
+      }
+    });
 
-    const std::shared_ptr<const MeshType> old_mesh0 = std::dynamic_pointer_cast<const MeshType>(old_mesh);
+  old_q = new_q;
+}
 
-    const std::shared_ptr<const MeshType> new_mesh0 = std::dynamic_pointer_cast<const MeshType>(new_mesh);
+template <size_t Dimension>
+void
+FluxingAdvectionSolver<Dimension>::_remapAllQuantities()
+{
+  const auto cell_to_face_matrix              = m_new_mesh->connectivity().cellToFaceMatrix();
+  const auto face_local_number_in_their_cells = m_new_mesh->connectivity().faceLocalNumbersInTheirCells();
 
-    FluxingAdvectionSolver<Dimension> solver(old_mesh0, new_mesh0);
+  MeshData<Dimension>& old_mesh_data = MeshDataManager::instance().getMeshData(*m_old_mesh);
 
-    FaceValue<double> fluxing_volumes(new_mesh0->connectivity());
-    fluxing_volumes.fill(0);
+  const CellValue<const double> old_Vj = old_mesh_data.Vj();
+  const CellValue<double> step_Vj      = copy(old_Vj);
 
-    std::vector<std::shared_ptr<const DiscreteFunctionVariant>> new_variables;
+  for (size_t jstep = 0; jstep < m_number_of_cycles; ++jstep) {
+    for (auto& remapped_q : m_remapped_list) {
+      std::visit([&](auto&& old_q) { this->_remapOne(step_Vj, old_q); }, remapped_q);
+    }
 
-    // for (auto&& variable_v : remapped_variables) {
-    //   std::visit(
-    //     [&](auto&& variable) {
-    //       using DiscreteFunctionT = std::decay_t<decltype(variable)>;
-    //       if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
-    //         remapUsingFluxing(new_mesh0, fluxing_volumes, variable);
-    //       }
-    //     },
-    //     variable_v->discreteFunction());
-    // }
+    parallel_for(
+      m_new_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        const auto& cell_to_face = cell_to_face_matrix[cell_id];
+        for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+          const FaceId face_id = cell_to_face[i_face];
+          CellId donnor_id     = m_donnor_cell[face_id];
 
-    return remapped_variables;   // std::make_shared<std::vector<std::shared_ptr<const
+          double flux = ((donnor_id == cell_id) ? -1 : 1) * m_cycle_fluxing_volume[face_id];
+          step_Vj[cell_id] += flux;
+        }
+      });
 
-    // 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");
-  }
+    CellValue<double> inv_Vj(m_old_mesh->connectivity());
+    parallel_for(
+      m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { inv_Vj[cell_id] = 1 / step_Vj[cell_id]; });
+
+    for (auto& remapped_q : m_remapped_list) {
+      std::visit(
+        [&](auto&& new_q) {
+          using CellDataType = std::decay_t<decltype(new_q)>;
+          static_assert(is_item_value_v<CellDataType> or is_item_array_v<CellDataType>, "invalid data type");
+
+          if constexpr (is_item_value_v<CellDataType>) {
+            parallel_for(
+              m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { new_q[cell_id] *= inv_Vj[cell_id]; });
+          } else if constexpr (is_item_array_v<CellDataType>) {
+            parallel_for(
+              m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+                auto array              = new_q[cell_id];
+                const double inv_volume = inv_Vj[cell_id];
+
+                for (size_t i = 0; i < array.size(); ++i) {
+                  array[i] *= inv_volume;
+                }
+              });
+          }
+        },
+        remapped_q);
+    }
   }
 }
 
-std::shared_ptr<const DiscreteFunctionVariant>
-FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
-                              const std::shared_ptr<const DiscreteFunctionVariant>& remapped_variable)
+template <size_t Dimension>
+std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
+FluxingAdvectionSolver<Dimension>::remap(
+  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& quantity_list)
 {
-  const std::shared_ptr<const IMesh> old_mesh = getCommonMesh({remapped_variable});
-
-  if (not checkDiscretizationType({remapped_variable}, DiscreteFunctionType::P0)) {
-    throw NormalError("acoustic solver expects P0 functions");
+  for (auto&& variable_v : quantity_list) {
+    std::visit(
+      [&](auto&& variable) {
+        using DiscreteFunctionT = std::decay_t<decltype(variable)>;
+        if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
+          this->_storeValues(variable);
+        } else {
+          throw UnexpectedError("incompatible mesh types");
+        }
+      },
+      variable_v->discreteFunction());
   }
 
-  switch (old_mesh->dimension()) {
-  case 1: {
-    throw NormalError("Not yet implemented in 1d");
-  }
-  case 2: {
-    constexpr size_t Dimension = 2;
-    using MeshType             = Mesh<Connectivity<Dimension>>;
+  this->_remapAllQuantities();
 
-    const std::shared_ptr<const MeshType> old_mesh0 = std::dynamic_pointer_cast<const MeshType>(old_mesh);
+  std::vector<std::shared_ptr<const DiscreteFunctionVariant>> new_variables;
 
-    const std::shared_ptr<const MeshType> new_mesh0 = std::dynamic_pointer_cast<const MeshType>(new_mesh);
-
-    FluxingAdvectionSolver<Dimension> solver(old_mesh0, new_mesh0);
-
-    FaceValue<double> fluxing_volumes = solver.computeFluxVolume();
-    const size_t number_of_cycles     = calculateRemapCycles(old_mesh0, fluxing_volumes);
-
-    std::cout << " number_of_cycle " << number_of_cycles << "\n";
-
-    DiscreteFunctionVariant new_variable = std::visit(
-      [&](auto&& variable) -> DiscreteFunctionVariant {
+  for (size_t i = 0; i < quantity_list.size(); ++i) {
+    std::visit(
+      [&](auto&& variable) {
         using DiscreteFunctionT = std::decay_t<decltype(variable)>;
+        using DataType          = std::decay_t<typename DiscreteFunctionT::data_type>;
+
         if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
-          return remapUsingFluxing(old_mesh0, new_mesh0, fluxing_volumes, number_of_cycles, variable);
+          if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+            new_variables.push_back(std::make_shared<DiscreteFunctionVariant>(
+              DiscreteFunctionT(m_new_mesh, std::get<CellValue<DataType>>(m_remapped_list[i]))));
+          } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+            new_variables.push_back(std::make_shared<DiscreteFunctionVariant>(
+              DiscreteFunctionT(m_new_mesh, std::get<CellArray<DataType>>(m_remapped_list[i]))));
+          } else {
+            throw UnexpectedError("invalid discrete function type");
+          }
         } else {
           throw UnexpectedError("incompatible mesh types");
         }
       },
-      remapped_variable->discreteFunction());
+      quantity_list[i]->discreteFunction());
+  }
+
+  return new_variables;
+}
 
-    return std::make_shared<DiscreteFunctionVariant>(new_variable);
+std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
+advectByFluxing(const std::shared_ptr<const IMesh> i_new_mesh,
+                const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& remapped_variables)
+{
+  if (not hasSameMesh(remapped_variables)) {
+    throw NormalError("remapped quantities are not defined on the same mesh");
+  }
+
+  const std::shared_ptr<const IMesh> i_old_mesh = getCommonMesh(remapped_variables);
+
+  switch (i_old_mesh->dimension()) {
+  case 1: {
+    return FluxingAdvectionSolver<1>{i_old_mesh, i_new_mesh}.remap(remapped_variables);
+  }
+  case 2: {
+    return FluxingAdvectionSolver<2>{i_old_mesh, i_new_mesh}.remap(remapped_variables);
   }
   case 3: {
-    throw NotImplementedError("Fluxing advection solver not implemented in dimension 3");
+    return FluxingAdvectionSolver<3>{i_old_mesh, i_new_mesh}.remap(remapped_variables);
   }
   default: {
     throw UnexpectedError("Invalid mesh dimension");
diff --git a/src/scheme/FluxingAdvectionSolver.hpp b/src/scheme/FluxingAdvectionSolver.hpp
index 43a3e6424..05dce5055 100644
--- a/src/scheme/FluxingAdvectionSolver.hpp
+++ b/src/scheme/FluxingAdvectionSolver.hpp
@@ -6,12 +6,8 @@
 
 #include <vector>
 
-std::vector<std::shared_ptr<const DiscreteFunctionVariant>> FluxingAdvectionSolverHandler(
+std::vector<std::shared_ptr<const DiscreteFunctionVariant>> advectByFluxing(
   const std::shared_ptr<const IMesh> new_mesh,
   const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& remapped_variables);
 
-std::shared_ptr<const DiscreteFunctionVariant> FluxingAdvectionSolverHandler(
-  const std::shared_ptr<const IMesh> new_mesh,
-  const std::shared_ptr<const DiscreteFunctionVariant>& remapped_variables);
-
 #endif   // FLUXING_ADVECION_SOLVER_HPP
-- 
GitLab