diff --git a/src/scheme/CellbyCellLimitation.hpp b/src/scheme/CellbyCellLimitation.hpp
index f9b8fc6028db2775fdaf09b74d1de9f5f04d5406..97d11a07e035595735861a173df6210179b4f2d8 100644
--- a/src/scheme/CellbyCellLimitation.hpp
+++ b/src/scheme/CellbyCellLimitation.hpp
@@ -17,12 +17,15 @@
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshTraits.hpp>
#include <mesh/MeshVariant.hpp>
+#include <mesh/StencilManager.hpp>
#include <scheme/DiscreteFunctionDPk.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionDPkVector.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
#include <scheme/PolynomialReconstruction.hpp>
+#include <scheme/PolynomialReconstructionDescriptor.hpp>
#include <utils/PugsTraits.hpp>
#include <variant>
@@ -39,6 +42,7 @@ class CellByCellLimitation
public:
void
density_limiter(const MeshType& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
const DiscreteFunctionP0<const double>& rho,
DiscreteFunctionDPk<Dimension, double>& rho_L) const
{
@@ -46,16 +50,20 @@ class CellByCellLimitation
const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
const auto& xr = mesh.xr();
- const auto& xl = mesh.xl();
+ // const auto& xl = mesh.xl();
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- auto stencil = StencilManager::instance().getStencil(mesh.connectivity(), 1);
-
- auto xj = mesh_data.xj();
+ // auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1);
+ auto stencil = StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
+ symmetry_boundary_descriptor_list);
+ auto xj = mesh_data.xj();
+ const auto& xl = mesh_data.xl();
- const QuadratureFormula<1> qf =
- QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
+ // const QuadratureFormula<1> qf =
+ // QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
@@ -81,21 +89,31 @@ class CellByCellLimitation
rho_bar_max = std::max(rho_bar_max, rho_xk);
}
- auto face_list = cell_to_face_matrix[cell_id];
- for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
- const FaceId face_id = face_list[i_face];
+ if constexpr (Dimension == 2) {
+ auto face_list = cell_to_face_matrix[cell_id];
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
- const Rd& x0 = xr[face_to_node_matrix[face_id][0]];
- const Rd& x1 = xl[face_id][0];
- const Rd& x2 = xr[face_to_node_matrix[face_id][1]];
+ const Rd& x0 = xr[face_to_node_matrix[face_id][0]];
+ const Rd& x2 = xr[face_to_node_matrix[face_id][1]];
+ const Rd& x1 = .5 * (x0 + x2);
- const LineParabolicTransformation<Dimension> t(x0, x1, x2);
- for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) {
- const double rho_xk = rho_L[cell_id](t(qf.point(iq)));
+ const double rho_x0 = rho_L[cell_id](x0);
+ const double rho_x1 = rho_L[cell_id](x1);
+ const double rho_x2 = rho_L[cell_id](x2);
- rho_bar_min = std::min(rho_bar_min, rho_xk);
- rho_bar_max = std::max(rho_bar_max, rho_xk);
+ rho_bar_min = std::min(rho_bar_min, std::min(rho_x0, std::min(rho_x1, rho_x2)));
+ rho_bar_max = std::max(rho_bar_max, std::max(rho_x0, std::max(rho_x1, rho_x2)));
}
+ // const LineParabolicTransformation<Dimension> t(x0, x1, x2);
+ // for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) {
+ // const double rho_xk = rho_L[cell_id](t(qf.point(iq)));
+
+ // rho_bar_min = std::min(rho_bar_min, rho_xk);
+ // rho_bar_max = std::max(rho_bar_max, rho_xk);
+ //}
+ } else {
}
const double eps = 1E-14;
@@ -122,11 +140,17 @@ class CellByCellLimitation
}
void
- specific_internal_nrj_limiter(const MeshType& mesh,
- const DiscreteFunctionP0<const double>& rho,
- const DiscreteFunctionDPk<Dimension, double>& rho_L
- DiscreteFunctionP0<const double>& epsilon,
- DiscreteFunctionDPk<Dimension, double>& epsilon_R) const
+ specific_internal_nrj_limiter(
+ const MeshType& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
+
+ // const DiscreteFunctionP0<const double>& rho,
+ // const DiscreteFunctionDPk<Dimension, double>& rho_L,
+ const DiscreteFunctionP0<const double>& epsilon,
+ // const DiscreteFunctionDPk<Dimension, double>&
+ auto epsilon_R(const CellId cell_id, const Rd& x),
+ CellValue<const double>& lambda_epsilon) // const
+ // DiscreteFunctionDPk<Dimension, double>& epsilon_R) const
{
const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
@@ -136,12 +160,15 @@ class CellByCellLimitation
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- auto stencil = StencilManager::instance().getStencil(mesh.connectivity(), 1);
-
+ // auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1);
+ auto stencil = StencilManager::instance()
+ .getCellToCellStencilArray(mesh.connectivity(),
+ StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
+ symmetry_boundary_descriptor_list);
auto xj = mesh_data.xj();
- const QuadratureFormula<1> qf =
- QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
+ // const QuadratureFormula<1> qf =
+ // QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
@@ -171,17 +198,23 @@ class CellByCellLimitation
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- const Rd& x0 = xr[face_to_node_matrix[face_id][0]];
- const Rd& x1 = xl[face_id][0];
- const Rd& x2 = xr[face_to_node_matrix[face_id][1]];
+ const Rd& x0 = xr[face_to_node_matrix[face_id][0]];
+ const Rd& x1 = xl[face_id][0];
+ const Rd& x2 = xr[face_to_node_matrix[face_id][1]];
+ const double epsilon_x0 = epsilon_R(cell_id, x0);
+ const double epsilon_x1 = epsilon_R(cell_id, x1);
+ const double epsilon_x2 = epsilon_R(cell_id, x2);
- const LineParabolicTransformation<Dimension> t(x0, x1, x2);
- for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) {
- const double epsilon_xk = epsilon_R(cell_id, t(qf.point(iq)));
+ epsilon_R_min = std::min(epsilon_R_min, std::min(epsilon_x0, std::min(epsilon_x1, epsilon_x2)));
+ epsilon_R_max = std::max(epsilon_R_max, std::max(epsilon_x0, std::max(epsilon_x1, epsilon_x2)));
- epsilon_R_min = std::min(epsilon_R_min, epsilon_xk);
- epsilon_R_max = std::max(epsilon_R_max, epsilon_xk);
- }
+ // const LineParabolicTransformation<Dimension> t(x0, x1, x2);
+ // for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) {
+ // const double epsilon_xk = epsilon_R(cell_id, t(qf.point(iq)));
+
+ // epsilon_R_min = std::min(epsilon_R_min, epsilon_xk);
+ // epsilon_R_max = std::max(epsilon_R_max, epsilon_xk);
+ //}
}
const double eps = 1E-14;