Skip to content
Snippets Groups Projects
Commit 71f55967 authored by Emmanuel Labourasse's avatar Emmanuel Labourasse
Browse files

Compute the fluxing volumes

parent e5ca1c71
No related branches found
No related tags found
1 merge request!167Improve fluxing based remapping
......@@ -590,13 +590,23 @@ SchemeModule::SchemeModule()
}
));
this->_addBuiltinFunction("fluxing_advection", std::function(
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>& old_q)
const std::shared_ptr<const DiscreteFunctionVariant>& remapped_variables)
-> std::shared_ptr<const DiscreteFunctionVariant> {
return FluxingAdvectionSolverHandler(new_mesh, old_q);
}));
return FluxingAdvectionSolverHandler(new_mesh, remapped_variables);
}
));
MathFunctionRegisterForVh{*this};
}
......
......@@ -4,10 +4,13 @@
#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
{
......@@ -15,10 +18,10 @@ class FluxingAdvectionSolver
using Rd = TinyVector<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()
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}
{}
~FluxingAdvectionSolver() = default;
};
template <>
FaceValue<double>
FluxingAdvectionSolver<1>::computeFluxVolume() const
{
throw NotImplementedError("Viens");
}
template <>
FaceValue<double>
FluxingAdvectionSolver<2>::computeFluxVolume() const
{
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");
}
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;
}
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}
{}
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();
~FluxingAdvectionSolver() = default;
};
DiscreteFunctionP0Vector<Dimension, double> new_q(new_mesh, copy(old_q.cellArrays()));
std::shared_ptr<const DiscreteFunctionVariant>
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);
return std::make_shared<DiscreteFunctionVariant>(solver.apply());
FaceValue<double> fluxing_volumes(old_mesh0->connectivity());
fluxing_volumes.fill(0);
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");
......
......@@ -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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment