diff --git a/src/scheme/CellbyCellLimitation.hpp b/src/scheme/CellbyCellLimitation.hpp
index 73a8ae06bec59c9b637c032bbb5f165c131fe15b..f9b8fc6028db2775fdaf09b74d1de9f5f04d5406 100644
--- a/src/scheme/CellbyCellLimitation.hpp
+++ b/src/scheme/CellbyCellLimitation.hpp
@@ -37,6 +37,168 @@ class CellByCellLimitation
using Rd = TinyVector<Dimension>;
public:
+ void
+ density_limiter(const MeshType& mesh,
+ const DiscreteFunctionP0<const double>& rho,
+ DiscreteFunctionDPk<Dimension, double>& rho_L) const
+ {
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+ const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+
+ const auto& xr = mesh.xr();
+ 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();
+
+ const QuadratureFormula<1> qf =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double rhoj = rho[cell_id];
+
+ double rho_min = rhoj;
+ double rho_max = rhoj;
+
+ const auto cell_stencil = stencil[cell_id];
+ for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ rho_min = std::min(rho_min, rho[cell_stencil[i_cell]]);
+ rho_max = std::max(rho_max, rho[cell_stencil[i_cell]]);
+ }
+
+ double rho_bar_min = rhoj;
+ double rho_bar_max = rhoj;
+
+ for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ const CellId cell_k_id = cell_stencil[i_cell];
+ const double rho_xk = rho_L[cell_id](xj[cell_k_id]);
+
+ rho_bar_min = std::min(rho_bar_min, rho_xk);
+ 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];
+
+ 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 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);
+ }
+ }
+
+ const double eps = 1E-14;
+ double coef1 = 1;
+ if (std::abs(rho_bar_max - rhoj) > eps) {
+ coef1 = (rho_max - rhoj) / ((rho_bar_max - rhoj));
+ }
+
+ double coef2 = 1.;
+ if (std::abs(rho_bar_min - rhoj) > eps) {
+ coef2 = (rho_min - rhoj) / ((rho_bar_min - rhoj));
+ }
+
+ const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+ auto coefficients = rho_L.coefficients(cell_id);
+
+ coefficients[0] = (1 - lambda) * rho[cell_id] + lambda * coefficients[0];
+
+ for (size_t i = 1; i < coefficients.size(); ++i) {
+ coefficients[i] *= lambda;
+ }
+ });
+ }
+
+ 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
+ {
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+ const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+
+ const auto& xr = mesh.xr();
+ 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();
+
+ const QuadratureFormula<1> qf =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double epsilonj = epsilon[cell_id];
+
+ double epsilon_min = epsilonj;
+ double epsilon_max = epsilonj;
+
+ const auto cell_stencil = stencil[cell_id];
+ for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ epsilon_min = std::min(epsilon_min, epsilon[cell_stencil[i_cell]]);
+ epsilon_max = std::max(epsilon_max, epsilon[cell_stencil[i_cell]]);
+ }
+
+ double epsilon_R_min = epsilonj;
+ double epsilon_R_max = epsilonj;
+
+ for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ const CellId cell_k_id = cell_stencil[i_cell];
+ const double epsilon_xk = epsilon_R(cell_id, xj[cell_k_id]);
+
+ epsilon_R_min = std::min(epsilon_R_min, epsilon_xk);
+ epsilon_R_max = std::max(epsilon_R_max, epsilon_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];
+
+ 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 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;
+ double coef1 = 1;
+ if (std::abs(epsilon_R_max - epsilonj) > eps) {
+ coef1 = (epsilon_max - epsilonj) / ((epsilon_R_max - epsilonj));
+ }
+
+ double coef2 = 1.;
+ if (std::abs(epsilon_R_min - epsilonj) > eps) {
+ coef2 = (epsilon_min - epsilonj) / ((epsilon_R_min - epsilonj));
+ }
+
+ lambda_epsilon[cell_id] = std::max(0., std::min(1., std::min(coef1, coef2)));
+ });
+ }
+
void
computeLimitorVolumicScalarQuantityMinModDukowiczGradient(const MeshType& mesh,
const CellValue<double>& q,