diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 9fc6d68bb0ae5295925cb6bcefce8308cbbf5cf2..4e6e873478e84299787aef9b4a2644af8406093b 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -590,13 +590,23 @@ SchemeModule::SchemeModule()
}
));
- this->_addBuiltinFunction("fluxing_advection", std::function(
-
- [](const std::shared_ptr<const IMesh> new_mesh,
- const std::shared_ptr<const DiscreteFunctionVariant>& old_q)
- -> std::shared_ptr<const DiscreteFunctionVariant> {
- return FluxingAdvectionSolverHandler(new_mesh, old_q);
- }));
+ 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);
+ }
+
+ ));
MathFunctionRegisterForVh{*this};
}
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index d3fcd4bb165e4c2b02c88157f97234e605eaf4a6..7ac8e85f21026ff9c6f83b26d5b949d7590972e8 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -4,21 +4,24 @@
#include <mesh/Connectivity.hpp>
#include <mesh/IMesh.hpp>
#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/SubItemValuePerItem.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 MeshType = Mesh<Connectivity<Dimension>>;
+ using MeshDataType = MeshData<Dimension>;
const std::shared_ptr<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>
@@ -72,32 +75,90 @@ class FluxingAdvectionSolver
// return Fnp1;
// }
- DiscreteFunctionP0<Dimension, double>
- apply()
- {
- 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;
- }
+ FaceValue<double> computeFluxVolume() const;
- FluxingAdvectionSolver(const std::shared_ptr<const MeshType> old_mesh,
- const std::shared_ptr<const MeshType> new_mesh,
- const DiscreteFunctionP0<Dimension, const double>& old_q)
+ 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}
{}
~FluxingAdvectionSolver() = default;
};
-std::shared_ptr<const DiscreteFunctionVariant>
+template <>
+FaceValue<double>
+FluxingAdvectionSolver<1>::computeFluxVolume() const
+{
+ throw NotImplementedError("Viens");
+}
+
+template <>
+FaceValue<double>
+FluxingAdvectionSolver<2>::computeFluxVolume() const
+{
+ if (m_new_mesh->shared_connectivity() != m_old_mesh->shared_connectivity()) {
+ throw NormalError("Old and new meshes must share the same connectivity");
+ }
+ 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());
+ 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);
+ });
+ return flux_volume;
+}
+
+template <>
+FaceValue<double>
+FluxingAdvectionSolver<3>::computeFluxVolume() const
+{
+ throw NotImplementedError("ViensViensViens");
+}
+
+template <typename MeshType, typename DataType>
+auto
+remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
+ [[maybe_unused]] const FaceValue<double>& fluxing_volumes,
+ const DiscreteFunctionP0<MeshType::Dimension, const DataType>& old_q)
+{
+ constexpr size_t Dimension = MeshType::Dimension;
+ // const Connectivity<Dimension>& connectivity = new_mesh->connectivity();
+
+ DiscreteFunctionP0<Dimension, DataType> new_q(new_mesh, copy(old_q.cellValues()));
+
+ return new_q;
+}
+
+template <typename MeshType>
+auto
+remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
+ [[maybe_unused]] const FaceValue<double>& fluxing_volumes,
+ const DiscreteFunctionP0Vector<MeshType::Dimension, const double>& old_q)
+{
+ constexpr size_t Dimension = MeshType::Dimension;
+ // const Connectivity<Dimension>& connectivity = new_mesh->connectivity();
+
+ DiscreteFunctionP0Vector<Dimension, double> new_q(new_mesh, copy(old_q.cellArrays()));
+
+ throw NotImplementedError("DiscreteFunctionP0Vector");
+
+ return new_q;
+}
+
+std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
- const std::shared_ptr<const DiscreteFunctionVariant>& old_q_v)
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& remapped_variables)
{
- const std::shared_ptr<const IMesh> old_mesh = getCommonMesh({old_q_v});
+ const std::shared_ptr<const IMesh> old_mesh = getCommonMesh(remapped_variables);
- if (not checkDiscretizationType({old_q_v}, DiscreteFunctionType::P0)) {
+ if (not checkDiscretizationType({remapped_variables}, DiscreteFunctionType::P0)) {
throw NormalError("acoustic solver expects P0 functions");
}
@@ -108,17 +169,116 @@ FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
const std::shared_ptr<const MeshType> old_mesh0 = std::dynamic_pointer_cast<const MeshType>(old_mesh);
- const DiscreteFunctionP0<Dimension, const double>& old_q =
- old_q_v->get<DiscreteFunctionP0<Dimension, const double>>();
+ 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(new_mesh0->connectivity());
+ fluxing_volumes.fill(0);
+
+ std::vector<std::shared_ptr<const DiscreteFunctionVariant>> new_variables;
+
+ // 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());
+ // }
+
+ return remapped_variables; // std::make_shared<std::vector<std::shared_ptr<const
+ // DiscreteFunctionVariant>>>(new_variables);
+ }
+ case 2: {
+ constexpr size_t Dimension = 2;
+
+ using MeshType = Mesh<Connectivity<Dimension>>;
+
+ const std::shared_ptr<const MeshType> old_mesh0 = std::dynamic_pointer_cast<const MeshType>(old_mesh);
+
+ 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(new_mesh0->connectivity());
+ fluxing_volumes.fill(0);
+
+ std::vector<std::shared_ptr<const DiscreteFunctionVariant>> new_variables;
+
+ // 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());
+ // }
+
+ return remapped_variables; // std::make_shared<std::vector<std::shared_ptr<const
+
+ // 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");
+ }
+ }
+}
+
+std::shared_ptr<const DiscreteFunctionVariant>
+FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
+ const std::shared_ptr<const DiscreteFunctionVariant>& remapped_variables)
+{
+ const std::shared_ptr<const IMesh> old_mesh = getCommonMesh({remapped_variables});
+
+ if (not checkDiscretizationType({remapped_variables}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ 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 std::shared_ptr<const MeshType> new_mesh0 = std::dynamic_pointer_cast<const MeshType>(new_mesh);
- FluxingAdvectionSolver<Dimension> solver(old_mesh0, new_mesh0, old_q);
+ // FluxingAdvectionSolver<Dimension> solver(old_mesh0, new_mesh0);
+
+ FaceValue<double> fluxing_volumes(old_mesh0->connectivity());
+ fluxing_volumes.fill(0);
- return std::make_shared<DiscreteFunctionVariant>(solver.apply());
+ return std::make_shared<DiscreteFunctionVariant>(
+ remapUsingFluxing(new_mesh0, fluxing_volumes,
+ remapped_variables->get<DiscreteFunctionP0<Dimension, const double>>()));
}
case 2: {
- throw NotImplementedError("Fluxing advection solver not implemented in dimension 2");
+ constexpr size_t Dimension = 2;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+
+ const std::shared_ptr<const MeshType> old_mesh0 = std::dynamic_pointer_cast<const MeshType>(old_mesh);
+
+ 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();
+ //(new_mesh0->connectivity());
+ // fluxing_volumes.fill(0);
+
+ return std::make_shared<DiscreteFunctionVariant>(
+ remapUsingFluxing(new_mesh0, fluxing_volumes,
+ remapped_variables->get<DiscreteFunctionP0<Dimension, const double>>()));
+
+ // throw NotImplementedError("Fluxing advection solver not implemented in dimension 2");
}
case 3: {
throw NotImplementedError("Fluxing advection solver not implemented in dimension 3");
diff --git a/src/scheme/FluxingAdvectionSolver.hpp b/src/scheme/FluxingAdvectionSolver.hpp
index d51ebb303a652d544419a1904c50785bf025c14c..43a3e6424432d6fc79a7cfcdbe9907a2130a8ee1 100644
--- a/src/scheme/FluxingAdvectionSolver.hpp
+++ b/src/scheme/FluxingAdvectionSolver.hpp
@@ -4,8 +4,14 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
+#include <vector>
+
+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);
+
std::shared_ptr<const DiscreteFunctionVariant> FluxingAdvectionSolverHandler(
const std::shared_ptr<const IMesh> new_mesh,
- const std::shared_ptr<const DiscreteFunctionVariant>& old_q);
+ const std::shared_ptr<const DiscreteFunctionVariant>& remapped_variables);
#endif // FLUXING_ADVECION_SOLVER_HPP