From 063dc78a1de1527b56044ee135e0b90acb0bba2d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Sat, 20 May 2023 10:03:51 +0200
Subject: [PATCH] Add boundary condition support for fluxing advection
Allowed bc tyes are
- inflow
- outflow
- symmetry
Inflow is not yet implemented for variables of type P0Vector
---
src/language/modules/SchemeModule.cpp | 17 +-
src/scheme/FluxingAdvectionSolver.cpp | 358 +++++++++++++++++++++++---
src/scheme/FluxingAdvectionSolver.hpp | 5 +
3 files changed, 343 insertions(+), 37 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 32628a2eb..813bde93a 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -677,16 +677,17 @@ 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 advectByFluxing(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 VariableBCDescriptor>>&
+ remapped_quantity_with_bc)
+ -> std::vector<std::shared_ptr<const DiscreteFunctionVariant>> {
+ return advectByFluxing(new_mesh, remapped_quantity_with_bc);
+ }
+
+ ));
MathFunctionRegisterForVh{*this};
}
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index 57f0d4842..97aad4bb8 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -1,6 +1,7 @@
#include <scheme/FluxingAdvectionSolver.hpp>
#include <language/utils/EvaluateAtPoints.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/IMesh.hpp>
#include <mesh/ItemArrayUtils.hpp>
@@ -9,10 +10,18 @@
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatFaceBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/SubItemValuePerItem.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/InflowBoundaryConditionDescriptor.hpp>
+#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+
+#include <variant>
+#include <vector>
template <size_t Dimension>
class FluxingAdvectionSolver
@@ -36,7 +45,25 @@ class FluxingAdvectionSolver
CellArray<double>>;
+ template <typename DataType>
+ class InflowValueBoundaryCondition;
+ class InflowArrayBoundaryCondition;
+ class OutflowBoundaryCondition;
+ class SymmetryBoundaryCondition;
+
+ using BoundaryConditionVariant = std::variant<InflowValueBoundaryCondition<double>, //
+ InflowValueBoundaryCondition<TinyVector<1>>,
+ InflowValueBoundaryCondition<TinyVector<2>>,
+ InflowValueBoundaryCondition<TinyVector<3>>,
+ InflowValueBoundaryCondition<TinyMatrix<1>>,
+ InflowValueBoundaryCondition<TinyMatrix<2>>,
+ InflowValueBoundaryCondition<TinyMatrix<3>>,
+ InflowArrayBoundaryCondition,
+ OutflowBoundaryCondition,
+ SymmetryBoundaryCondition>;
+
std::vector<RemapVariant> m_remapped_list;
+ std::vector<std::vector<BoundaryConditionVariant>> m_boundary_condition_list;
FaceValue<const CellId> m_donnor_cell;
FaceValue<const double> m_cycle_fluxing_volume;
@@ -55,21 +82,98 @@ class FluxingAdvectionSolver
m_remapped_list.emplace_back(copy(old_q.cellValues()));
}
- template <typename DataType>
void
- _storeValues(const DiscreteFunctionP0Vector<Dimension, const DataType>& old_q)
+ _storeValues(const DiscreteFunctionP0Vector<Dimension, const double>& old_q)
{
m_remapped_list.emplace_back(copy(old_q.cellArrays()));
}
+ template <typename DataType>
+ void
+ _storeValueBCList(const size_t i_quantity,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_list)
+ {
+ const Mesh<Connectivity<Dimension>>& mesh = *m_old_mesh;
+
+ for (auto& i_bc : bc_list) {
+ switch (i_bc->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ m_boundary_condition_list[i_quantity].emplace_back(
+ SymmetryBoundaryCondition(getMeshFlatFaceBoundary(mesh, i_bc->boundaryDescriptor())));
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::outflow: {
+ m_boundary_condition_list[i_quantity].emplace_back(
+ OutflowBoundaryCondition(getMeshFaceBoundary(mesh, i_bc->boundaryDescriptor())));
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::inflow: {
+ const InflowBoundaryConditionDescriptor& inflow_bc_descriptor =
+ dynamic_cast<const InflowBoundaryConditionDescriptor&>(*i_bc);
+
+ MeshFaceBoundary<Dimension> mesh_face_boundary =
+ getMeshFaceBoundary(mesh, inflow_bc_descriptor.boundaryDescriptor());
+
+ FaceValue<const Rd> xl = MeshDataManager::instance().getMeshData(*m_old_mesh).xl();
+ Array<const DataType> value_list =
+ InterpolateItemValue<DataType(Rd)>::template interpolate<ItemType::face>(inflow_bc_descriptor
+ .functionSymbolId(),
+ xl, mesh_face_boundary.faceList());
+
+ m_boundary_condition_list[i_quantity].emplace_back(
+ InflowValueBoundaryCondition<DataType>{mesh_face_boundary, value_list});
+ break;
+ }
+ default: {
+ std::ostringstream os;
+ os << "invalid boundary condition for advection: " << *i_bc;
+ throw NotImplementedError("inflow boundary condition");
+ }
+ }
+ }
+ }
+
+ void
+ _storeArrayBCList(const size_t i_quantity,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_list)
+ {
+ const Mesh<Connectivity<Dimension>>& mesh = *m_old_mesh;
+
+ for (auto& i_bc : bc_list) {
+ switch (i_bc->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ m_boundary_condition_list[i_quantity].emplace_back(
+ SymmetryBoundaryCondition(getMeshFlatFaceBoundary(mesh, i_bc->boundaryDescriptor())));
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::outflow: {
+ m_boundary_condition_list[i_quantity].emplace_back(
+ OutflowBoundaryCondition(getMeshFaceBoundary(mesh, i_bc->boundaryDescriptor())));
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::inflow: {
+ throw NotImplementedError("inflow boundary condition");
+ break;
+ }
+ default: {
+ std::ostringstream os;
+ os << "invalid boundary condition for advection: " << *i_bc;
+ throw NotImplementedError("inflow boundary condition");
+ }
+ }
+ }
+ }
+
template <typename CellDataType>
- void _remapOne(const CellValue<const double>& step_Vj, CellDataType& old_q);
+ void _remapOne(const CellValue<const double>& step_Vj,
+ const std::vector<BoundaryConditionVariant>& q_bc_list,
+ CellDataType& old_q);
void _remapAllQuantities();
public:
std::vector<std::shared_ptr<const DiscreteFunctionVariant>> //
- remap(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& quantities);
+ remap(const std::vector<std::shared_ptr<const VariableBCDescriptor>>& quantity_list_with_bc);
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)},
@@ -102,16 +206,16 @@ FluxingAdvectionSolver<Dimension>::_computeDonorCells(FaceValue<const double> al
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];
+ const auto& face_to_cell = face_to_cell_matrix[face_id];
+ const size_t i_face_in_cell = face_local_number_in_their_cells[face_id][0];
+ const CellId cell_id = face_to_cell[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 {
- 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 {
+ if (face_to_cell.size() == 2) {
donnor_cell[face_id] = face_to_cell[1];
+ } else {
+ donnor_cell[face_id] = std::numeric_limits<CellId::base_type>::max();
}
}
});
@@ -128,14 +232,20 @@ FluxingAdvectionSolver<1>::_computeDonorCells(FaceValue<const double> algebraic_
const auto face_to_cell_matrix = m_new_mesh->connectivity().faceToCellMatrix();
const auto cell_to_face_matrix = m_new_mesh->connectivity().cellToFaceMatrix();
+ 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];
+ const CellId cell_id = face_to_cell[0];
if (face_to_cell.size() == 1) {
- donnor_cell[face_id] = face_to_cell[0];
+ if ((algebraic_fluxing_volumes[face_id] <= 0) xor (cell_to_face_matrix[cell_id][1] == face_id)) {
+ donnor_cell[face_id] = std::numeric_limits<CellId::base_type>::max();
+ } else {
+ donnor_cell[face_id] = cell_id;
+ }
} 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 {
@@ -299,8 +409,9 @@ FluxingAdvectionSolver<Dimension>::_computeCycleNumber(FaceValue<double> fluxing
CellValue<size_t> ratio(m_old_mesh->connectivity());
parallel_for(
- m_old_mesh->numberOfCells(),
- PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = std::ceil(total_negative_flux[cell_id] / Vj[cell_id]); });
+ m_old_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ ratio[cell_id] = (1. / 0.6) * std::ceil(total_negative_flux[cell_id] / Vj[cell_id]);
+ });
synchronize(ratio);
size_t number_of_cycles = max(ratio);
@@ -329,8 +440,12 @@ FluxingAdvectionSolver<Dimension>::_computeGeometricalData()
template <size_t Dimension>
template <typename CellDataType>
void
-FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step_Vj, CellDataType& old_q)
+FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step_Vj,
+ const std::vector<BoundaryConditionVariant>& q_bc_list,
+ CellDataType& old_q)
{
+ using DataType = std::decay_t<typename CellDataType::data_type>;
+
static_assert(is_item_value_v<CellDataType> or is_item_array_v<CellDataType>, "invalid data type");
const auto cell_to_face_matrix = m_new_mesh->connectivity().cellToFaceMatrix();
@@ -353,6 +468,7 @@ FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step
});
}
+ // First we deal with inner faces
parallel_for(
m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
const auto& cell_to_face = cell_to_face_matrix[cell_id];
@@ -384,6 +500,97 @@ FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step
}
});
+ // Now we deal with boundary faces
+ for (const auto& bc_descriptor : q_bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using TypeOfBC = std::decay_t<decltype(bc)>;
+
+ if constexpr (std::is_same_v<TypeOfBC, SymmetryBoundaryCondition>) {
+ auto face_list = bc.faceList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const double fluxing_volume = m_cycle_fluxing_volume[face_id];
+ const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+ if (fluxing_volume > 1E-12 * step_Vj[face_cell_id]) {
+ std::ostringstream error_msg;
+ error_msg << "invalid symmetry for face " << face_id
+ << " (number=" << m_old_mesh->connectivity().faceNumber()[face_id] << ")\n"
+ << "face has non zero fluxing volume " << fluxing_volume
+ << " (cell volume = " << step_Vj[face_cell_id] << ")\n";
+ throw NormalError(error_msg.str());
+ }
+ }
+ } else if constexpr (std::is_same_v<TypeOfBC, OutflowBoundaryCondition>) {
+ auto face_list = bc.faceList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const CellId cell_id = m_donnor_cell[face_id];
+ const double fluxing_volume = m_cycle_fluxing_volume[face_id];
+ if (cell_id != std::numeric_limits<CellId::base_type>::max()) {
+ if constexpr (is_item_array_v<CellDataType>) {
+ for (size_t i = 0; i < new_q[cell_id].size(); ++i) {
+ new_q[cell_id][i] -= fluxing_volume * old_q[cell_id][i];
+ }
+ } else {
+ new_q[cell_id] -= fluxing_volume * old_q[cell_id];
+ }
+ } else {
+ const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+ if (fluxing_volume > 1E-12 * step_Vj[face_cell_id]) {
+ std::ostringstream error_msg;
+ error_msg << "invalid outflow for face " << face_id
+ << " (number=" << m_old_mesh->connectivity().faceNumber()[face_id] << ")"
+ << "face has inflow fluxing volume " << fluxing_volume
+ << " (cell volume = " << step_Vj[face_cell_id] << ")\n";
+ throw NormalError(error_msg.str());
+ }
+ }
+ }
+ } else if constexpr (std::is_same_v<TypeOfBC, InflowValueBoundaryCondition<DataType>>) {
+ if constexpr (is_item_value_v<CellDataType>) {
+ auto face_list = bc.faceList();
+ auto value_list = bc.valueList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const CellId cell_id = m_donnor_cell[face_id];
+ const double fluxing_volume = m_cycle_fluxing_volume[face_id];
+ if (cell_id == std::numeric_limits<CellId::base_type>::max()) {
+ if constexpr (is_item_value_v<CellDataType>) {
+ Assert(face_to_cell_matrix[face_id].size() == 1, "boundary face must be connected to a single cell");
+ const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+ new_q[face_cell_id] += fluxing_volume * value_list[i_face];
+ } else {
+ throw UnexpectedError("invalid boundary condition for fluxing advection");
+ }
+ } else {
+ const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+ if (fluxing_volume > 1E-12 * step_Vj[face_cell_id]) {
+ std::ostringstream error_msg;
+ error_msg << "invalid inflow for face " << face_id
+ << " (number=" << m_old_mesh->connectivity().faceNumber()[face_id] << ")"
+ << "face has outflow fluxing volume " << fluxing_volume
+ << " (cell volume = " << step_Vj[face_cell_id] << ")\n";
+ throw NormalError(error_msg.str());
+ }
+ }
+ }
+ } else {
+ throw UnexpectedError("invalid boundary condition for fluxing advection");
+ }
+ } else if constexpr (std::is_same_v<TypeOfBC, InflowArrayBoundaryCondition>) {
+ if constexpr (is_item_array_v<CellDataType>) {
+ throw NotImplementedError("advect array");
+ } else {
+ throw UnexpectedError("invalid boundary condition for fluxing advection");
+ }
+ } else {
+ throw UnexpectedError("invalid boundary condition for fluxing advection");
+ }
+ },
+ bc_descriptor);
+ }
+
synchronize(new_q);
old_q = new_q;
}
@@ -399,10 +606,12 @@ FluxingAdvectionSolver<Dimension>::_remapAllQuantities()
const CellValue<const double> old_Vj = old_mesh_data.Vj();
const CellValue<double> step_Vj = copy(old_Vj);
-
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 (size_t i = 0; i < m_remapped_list.size(); ++i) {
+ Assert(m_remapped_list.size() == m_boundary_condition_list.size());
+ auto& remapped_q = m_remapped_list[i];
+ auto& q_bc_list = m_boundary_condition_list[i];
+ std::visit([&](auto&& old_q) { this->_remapOne(step_Vj, q_bc_list, old_q); }, remapped_q);
}
parallel_for(
@@ -452,26 +661,32 @@ FluxingAdvectionSolver<Dimension>::_remapAllQuantities()
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::vector<std::shared_ptr<const VariableBCDescriptor>>& quantity_list_with_bc)
{
- for (auto&& variable_v : quantity_list) {
+ m_boundary_condition_list.resize(quantity_list_with_bc.size());
+ for (size_t i = 0; i < quantity_list_with_bc.size(); ++i) {
+ const auto& quantity_v_with_bc = quantity_list_with_bc[i];
+ const auto& quantity_v = quantity_v_with_bc->discreteFunctionVariant();
+ const auto& bc_list = quantity_v_with_bc->bcDescriptorList();
std::visit(
[&](auto&& variable) {
using DiscreteFunctionT = std::decay_t<decltype(variable)>;
if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
this->_storeValues(variable);
+ using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+ this->_storeValueBCList<DataType>(i, bc_list);
} else {
throw UnexpectedError("incompatible mesh types");
}
},
- variable_v->discreteFunction());
+ quantity_v->discreteFunction());
}
this->_remapAllQuantities();
std::vector<std::shared_ptr<const DiscreteFunctionVariant>> new_variables;
- for (size_t i = 0; i < quantity_list.size(); ++i) {
+ for (size_t i = 0; i < quantity_list_with_bc.size(); ++i) {
std::visit(
[&](auto&& variable) {
using DiscreteFunctionT = std::decay_t<decltype(variable)>;
@@ -491,16 +706,101 @@ FluxingAdvectionSolver<Dimension>::remap(
throw UnexpectedError("incompatible mesh types");
}
},
- quantity_list[i]->discreteFunction());
+ quantity_list_with_bc[i]->discreteFunctionVariant()->discreteFunction());
}
return new_variables;
}
+template <size_t Dimension>
+template <typename DataType>
+class FluxingAdvectionSolver<Dimension>::InflowValueBoundaryCondition
+{
+ private:
+ const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+
+ Array<const DataType> m_value_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ const Array<const DataType>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ InflowValueBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary, Array<const DataType> value_list)
+ : m_mesh_face_boundary(mesh_face_boundary), m_value_list(value_list)
+ {
+ ;
+ }
+
+ ~InflowValueBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class FluxingAdvectionSolver<Dimension>::InflowArrayBoundaryCondition
+{
+};
+
+template <size_t Dimension>
+class FluxingAdvectionSolver<Dimension>::OutflowBoundaryCondition
+{
+ private:
+ const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ OutflowBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary)
+ : m_mesh_face_boundary(mesh_face_boundary)
+ {
+ ;
+ }
+
+ ~OutflowBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class FluxingAdvectionSolver<Dimension>::SymmetryBoundaryCondition
+{
+ private:
+ const MeshFlatFaceBoundary<Dimension> m_mesh_flat_face_boundary;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_flat_face_boundary.faceList();
+ }
+
+ SymmetryBoundaryCondition(const MeshFlatFaceBoundary<Dimension>& mesh_flat_face_boundary)
+ : m_mesh_flat_face_boundary(mesh_flat_face_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
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)
+ const std::vector<std::shared_ptr<const VariableBCDescriptor>>& remapped_variables_with_bc)
{
+ std::vector<std::shared_ptr<const DiscreteFunctionVariant>> remapped_variables;
+ for (auto i_remapped_var_with_bc : remapped_variables_with_bc) {
+ remapped_variables.push_back(i_remapped_var_with_bc->discreteFunctionVariant());
+ }
+
if (not hasSameMesh(remapped_variables)) {
throw NormalError("remapped quantities are not defined on the same mesh");
}
@@ -509,13 +809,13 @@ advectByFluxing(const std::shared_ptr<const IMesh> i_new_mesh,
switch (i_old_mesh->dimension()) {
case 1: {
- return FluxingAdvectionSolver<1>{i_old_mesh, i_new_mesh}.remap(remapped_variables);
+ return FluxingAdvectionSolver<1>{i_old_mesh, i_new_mesh}.remap(remapped_variables_with_bc);
}
case 2: {
- return FluxingAdvectionSolver<2>{i_old_mesh, i_new_mesh}.remap(remapped_variables);
+ return FluxingAdvectionSolver<2>{i_old_mesh, i_new_mesh}.remap(remapped_variables_with_bc);
}
case 3: {
- return FluxingAdvectionSolver<3>{i_old_mesh, i_new_mesh}.remap(remapped_variables);
+ return FluxingAdvectionSolver<3>{i_old_mesh, i_new_mesh}.remap(remapped_variables_with_bc);
}
default: {
throw UnexpectedError("Invalid mesh dimension");
diff --git a/src/scheme/FluxingAdvectionSolver.hpp b/src/scheme/FluxingAdvectionSolver.hpp
index 05dce5055..2cbb653cb 100644
--- a/src/scheme/FluxingAdvectionSolver.hpp
+++ b/src/scheme/FluxingAdvectionSolver.hpp
@@ -3,6 +3,7 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/VariableBCDescriptor.hpp>
#include <vector>
@@ -10,4 +11,8 @@ 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::vector<std::shared_ptr<const DiscreteFunctionVariant>> advectByFluxing(
+ const std::shared_ptr<const IMesh> new_mesh,
+ const std::vector<std::shared_ptr<const VariableBCDescriptor>>& remapped_variables_with_bc);
+
#endif // FLUXING_ADVECION_SOLVER_HPP
--
GitLab