From f0f82533a02f81cedfa2bbd7f975c85aebe325f3 Mon Sep 17 00:00:00 2001
From: Alexandre Gangloff <alexandre.gangloff@cea.fr>
Date: Fri, 7 Mar 2025 08:37:31 +0100
Subject: [PATCH] Clean up the order 2 hyperelastic monolithic solver code

---
 src/scheme/LocalDtUtils.hpp             | 124 +++++++
 src/scheme/Order2HyperelasticSolver.cpp | 437 +-----------------------
 src/scheme/Order2Limiters.hpp           | 330 ++++++++++++++++++
 3 files changed, 460 insertions(+), 431 deletions(-)
 create mode 100644 src/scheme/LocalDtUtils.hpp
 create mode 100644 src/scheme/Order2Limiters.hpp

diff --git a/src/scheme/LocalDtUtils.hpp b/src/scheme/LocalDtUtils.hpp
new file mode 100644
index 000000000..cfd8ab7ff
--- /dev/null
+++ b/src/scheme/LocalDtUtils.hpp
@@ -0,0 +1,124 @@
+#ifndef LOCAL_DT_UTILS_HPP
+#define LOCAL_DT_UTILS_HPP
+
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/StencilArray.hpp>
+#include <mesh/StencilManager.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <scheme/PolynomialReconstructionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Socket.hpp>
+
+template <MeshConcept MeshType>
+std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
+             const std::shared_ptr<const DiscreteFunctionVariant>,
+             const std::shared_ptr<const DiscreteFunctionVariant>,
+             const std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_state_law(const size_t& law,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& CG_v,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& lambda_v,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& mu_v,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& gamma_v,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& p_inf_v)
+  {
+    std::shared_ptr mesh_v = getCommonMesh({rho_v, u_v});
+
+    const DiscreteFunctionP0<const double>& rho_d                         = rho_v->get<DiscreteFunctionP0<const double>>();
+    const DiscreteFunctionP0<const TinyVector<MeshType::Dimension>>& u_d  = u_v->get<DiscreteFunctionP0<const TinyVector<MeshType::Dimension>>>();
+    const DiscreteFunctionP0<const double>& E_d                           = E_v->get<DiscreteFunctionP0<const double>>();
+    const DiscreteFunctionP0<const TinyMatrix<MeshType::Dimension>>& CG_d = CG_v->get<DiscreteFunctionP0<const TinyMatrix<MeshType::Dimension>>>();
+    const DiscreteFunctionP0<const double>& lambda_d                      = lambda_v->get<DiscreteFunctionP0<const double>>();
+    const DiscreteFunctionP0<const double>& mu_d                          = mu_v->get<DiscreteFunctionP0<const double>>();
+    const DiscreteFunctionP0<const double>& gamma_d                       = gamma_v->get<DiscreteFunctionP0<const double>>();
+    const DiscreteFunctionP0<const double>& p_inf_d                       = p_inf_v->get<DiscreteFunctionP0<const double>>();
+
+    const TinyMatrix<MeshType::Dimension> I = identity;
+    const DiscreteFunctionP0<const double>& eps_d = E_d - 0.5 * dot(u_d, u_d);
+
+    DiscreteFunctionP0<double> p_d(mesh_v);
+    DiscreteFunctionP0<TinyMatrix<MeshType::Dimension>> sigma_d(mesh_v);
+    DiscreteFunctionP0<double> aL_d(mesh_v);
+    DiscreteFunctionP0<double> aT_d(mesh_v);
+
+    if (law == 1) {
+      p_d = (gamma_d - 1) * rho_d * eps_d - p_inf_d * gamma_d;
+      sigma_d = (mu_d / sqrt(det(CG_d)) * (CG_d - I) + lambda_d / sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I) - p_d * I;
+      aL_d = sqrt((lambda_d + 2 * mu_d) / rho_d + gamma_d * (p_d + p_inf_d) / rho_d);
+      aT_d = sqrt(mu_d / rho_d);
+    } else if (law == 2){
+      p_d = (gamma_d - 1) * rho_d * eps_d - gamma_d * p_inf_d;
+      sigma_d = mu_d / sqrt(det(CG_d)) * (CG_d - 1. / 3. * trace(CG_d) * I) - p_d * I;
+      aL_d = sqrt((2 * mu_d) / rho_d + (gamma_d * (p_d + p_inf_d)) / rho_d);
+      aT_d = sqrt(mu_d / rho_d);
+    } else {
+      sigma_d = (mu_d / sqrt(det(CG_d)) * (CG_d - I) + lambda_d / sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I);
+      p_d = -1./3. * trace(sigma_d);
+      aL_d = sqrt((lambda_d + 2 * mu_d) / rho_d);
+      aT_d = sqrt(mu_d / rho_d);
+    }
+    return {std::make_shared<DiscreteFunctionVariant>(sigma_d), 
+            std::make_shared<DiscreteFunctionVariant>(p_d),
+            std::make_shared<DiscreteFunctionVariant>(aL_d),
+            std::make_shared<DiscreteFunctionVariant>(aT_d)};
+  }
+
+  template<MeshConcept MeshType>
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  mesh_correction(const std::shared_ptr<const MeshVariant>& i_mesh1,
+                  const std::shared_ptr<const MeshVariant>& i_mesh2,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& CG,
+                  std::vector<NodeId>& map1,
+                  std::vector<NodeId>& map2)
+  {
+    const MeshType mesh1 = *i_mesh1->get<MeshType>();
+    const MeshType mesh2 = *i_mesh2->get<MeshType>();
+
+    NodeValue<TinyVector<MeshType::Dimension>> xr1     = copy(mesh1.xr());
+    NodeValue<TinyVector<MeshType::Dimension>> new_xr2 = copy(mesh2.xr());
+
+    size_t n = map1.size();
+
+    for (size_t i = 0; i < n; i++) {
+      for (size_t j = 0; j < MeshType::Dimension; j++) {
+        new_xr2[map2[i]][j] = xr1[map1[i]][j];
+      }
+    }
+
+    std::shared_ptr<const MeshType> new_mesh2 = std::make_shared<MeshType>(mesh2.shared_connectivity(), new_xr2);
+
+    const std::shared_ptr<const DiscreteFunctionVariant>& new_rho = shallowCopy(new_mesh2, rho);
+    const std::shared_ptr<const DiscreteFunctionVariant>& new_u   = shallowCopy(new_mesh2, u);
+    const std::shared_ptr<const DiscreteFunctionVariant>& new_E   = shallowCopy(new_mesh2, E);
+    const std::shared_ptr<const DiscreteFunctionVariant>& new_CG  = shallowCopy(new_mesh2, CG);
+
+    return std::make_tuple(new_mesh2,new_rho,new_u,new_E,new_CG);
+  }
+
+#endif //LOCAL_DT_UTILS_HPP
\ No newline at end of file
diff --git a/src/scheme/Order2HyperelasticSolver.cpp b/src/scheme/Order2HyperelasticSolver.cpp
index 97ef40f91..580119d6e 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -21,6 +21,8 @@
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/LocalDtUtils.hpp>
+#include <scheme/Order2Limiters.hpp>
 #include <scheme/PolynomialReconstruction.hpp>
 #include <scheme/PolynomialReconstructionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
@@ -521,433 +523,6 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
     return F;
   }
 
-void
-  _scalar_limiter(const MeshType& mesh,
-                  const DiscreteFunctionP0<const double>& f,
-                  DiscreteFunctionDPk<Dimension, double>& DPk_fh) const
-  {
-    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
-    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
-    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
-    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
-
-    const auto xj = mesh_data.xj();
-
-    CellValue<double> alph_J2{mesh.connectivity()};
-
-    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-      const double fj = f[cell_id];
-
-      double f_min = fj;
-      double f_max = fj;
-
-      const auto cell_stencil = stencil[cell_id];
-      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-        f_min = std::min(f_min, f[cell_stencil[i_cell]]);
-        f_max = std::max(f_max, f[cell_stencil[i_cell]]);
-      }
-
-      double f_bar_min = fj;
-      double f_bar_max = fj;
-
-      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-        const CellId cell_k_id = cell_stencil[i_cell];
-        const double f_xk    = DPk_fh[cell_id](xj[cell_k_id]);
-
-        f_bar_min = std::min(f_bar_min, f_xk);
-        f_bar_max = std::max(f_bar_max, f_xk);
-      }
-
-      const double eps = 1E-14;
-      double coef1     = 1;
-      if (std::abs(f_bar_max - fj) > eps) {
-        coef1 = (f_max - fj) / ((f_bar_max - fj));
-      }
-
-      double coef2 = 1.;
-      if (std::abs(f_bar_min - fj) > eps) {
-        coef2 = (f_min - fj) / ((f_bar_min - fj));
-      }
-
-      const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
-      alph_J2[cell_id] = lambda;
-
-      auto coefficients = DPk_fh.coefficients(cell_id);
-
-      coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
-      for (size_t i = 1; i < coefficients.size(); ++i) {
-        coefficients[i] *= lambda;
-      }
-    });
-  }
-
-  void
-  _vector_limiter(const MeshType& mesh,
-                  const DiscreteFunctionP0<const Rd>& f,
-                  DiscreteFunctionDPk<Dimension, Rd> DPk_fh,
-                  const CellValue<const Rdxd>& grad_u) const 
-  {
-    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
-    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
-    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
-    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
-                                                                        symmetry_boundary_descriptor_list);
-    
-    const auto xj = mesh_data.xj();
-
-    CellValue<Rdxd> sym_grad_u{mesh.connectivity()};
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId j){
-        sym_grad_u[j] = (grad_u[j] + transpose(grad_u[j]));
-    });
-
-    CellValue<SmallArray<double>> eigunvalues{mesh.connectivity()};
-    CellValue<std::vector<SmallVector<double>>> eigunvectors{mesh.connectivity()};
-    for(CellId j = 0; j < mesh.numberOfCells(); ++j){
-      EigenvalueSolver{}.computeForSymmetricMatrix(sym_grad_u[j],eigunvalues[j],eigunvectors[j]);
-    }
-
-    CellValue<Rd> n{mesh.connectivity()};
-    CellValue<Rd> t{mesh.connectivity()};
-    CellValue<Rd> l{mesh.connectivity()};
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-        n[cell_id] = zero;
-        t[cell_id] = zero;
-        l[cell_id] = zero;
-        //if(frobeniusNorm(sym_grad_u[cell_id]) > 1E-5){
-          if constexpr(Dimension == 2){
-            Rd nj;
-            Rd tj;
-            for(size_t i = 0; i < Dimension; ++i){
-              nj[i] = eigunvectors[cell_id][0][i];
-              tj[i] = eigunvectors[cell_id][1][i];
-            }
-            n[cell_id] = nj;
-            t[cell_id] = tj;
-            if(l2Norm(nj) != 0){
-              n[cell_id] = (1/l2Norm(nj))*n[cell_id];
-            }
-            if(l2Norm(tj) != 0){
-              t[cell_id] = (1/l2Norm(tj))*t[cell_id];
-            }
-          } else {
-            if constexpr(Dimension == 3){
-              Rd nj;
-              Rd tj;
-              Rd lj;
-              for(size_t i = 0; i < Dimension; ++i){
-                nj[i] = eigunvectors[cell_id][0][i];
-                tj[i] = eigunvectors[cell_id][1][i];
-                lj[i] = eigunvectors[cell_id][2][i];
-              }
-              n[cell_id] = nj;
-              t[cell_id] = tj;
-              l[cell_id] = lj;
-              if(l2Norm(nj) != 0){
-              n[cell_id] = (1/l2Norm(nj))*n[cell_id];
-              }
-              if(l2Norm(tj) != 0){
-                t[cell_id] = (1/l2Norm(tj))*t[cell_id];
-              }
-              if(l2Norm(lj) != 0){
-                l[cell_id] = (1/l2Norm(lj))*l[cell_id];
-              }
-            } else {
-                Rd nj;
-                nj[0] = 1.;
-                n[cell_id] = nj;
-            }
-          }
-        //}
-    });
-    
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-        const double fn = dot(f[cell_id],n[cell_id]);
-        const double ft = dot(f[cell_id],t[cell_id]);
-        const double fl = dot(f[cell_id],l[cell_id]);
-
-        double fn_min = fn;
-        double fn_max = fn;
-        double ft_min = ft;
-        double ft_max = ft;
-        double fl_min = fl;
-        double fl_max = fl;
-
-        const auto cell_stencil = stencil[cell_id];
-        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-          const double fn_k = dot(f[cell_stencil[i_cell]],n[cell_id]);
-          fn_min            = std::min(fn_min,fn_k);
-          fn_max            = std::max(fn_max,fn_k);
-
-          const double ft_k = dot(f[cell_stencil[i_cell]],t[cell_id]);
-          ft_min            = std::min(ft_min,ft_k);
-          ft_max            = std::max(ft_max,ft_k);
-
-          const double fl_k = dot(f[cell_stencil[i_cell]],l[cell_id]);
-          fl_min            = std::min(fl_min,fl_k);
-          fl_max            = std::max(fl_max,fl_k);
-        }
-
-        double fn_bar_min = fn;
-        double fn_bar_max = fn;
-        double ft_bar_min = ft;
-        double ft_bar_max = ft;
-        double fl_bar_min = fl;
-        double fl_bar_max = fl;
-
-        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-          const CellId cell_k_id = cell_stencil[i_cell];
-          const Rd f_xk          = DPk_fh[cell_id](xj[cell_k_id]);
-
-          const double fn_xk = dot(f_xk,n[cell_id]);
-          fn_bar_min         = std::min(fn_bar_min,fn_xk);
-          fn_bar_max         = std::max(fn_bar_max,fn_xk);
-
-          const double ft_xk = dot(f_xk,t[cell_id]);
-          ft_bar_min         = std::min(ft_bar_min,ft_xk);
-          ft_bar_max         = std::max(ft_bar_max,ft_xk);
-
-          const double fl_xk = dot(f_xk,l[cell_id]);
-          fl_bar_min         = std::min(fl_bar_min,fl_xk);
-          fl_bar_max         = std::max(fl_bar_max,fl_xk);
-        }
-
-        const double eps = 1E-14;
-        double coef1_n   = 1;
-        if(std::abs(fn_bar_max - fn) > eps){
-          coef1_n = (fn_max - fn) / (fn_bar_max - fn);
-        } 
-        double coef2_n = 1;
-        if(std::abs(fn_bar_min - fn) > eps){
-          coef2_n = (fn_min - fn) / (fn_bar_min - fn);
-        }
-        
-        double coef1_t = 1;
-        if(std::abs(ft_bar_max - ft) > eps){
-          coef1_t = (ft_max - ft) / (ft_bar_max - ft);
-        }
-        double coef2_t = 1;
-        if(std::abs(ft_bar_min - ft) > eps){
-          coef2_t = (ft_min - ft) / (ft_bar_min - ft);
-        }
-
-        double coef1_l = 1;
-        if(std::abs(fl_bar_max - fl) > eps){
-          coef1_l = (fl_max - fl) / (fl_bar_max - fl);
-        }
-        double coef2_l = 1;
-        if(std::abs(fl_bar_min - fl) > eps){
-          coef2_l = (fl_min - fl) / (fl_bar_min - fl);
-        }
-
-        const double lambda_n = std::max(0., std::min(1., std::min(coef1_n,coef2_n)));
-        const double lambda_t = std::max(0., std::min(1., std::min(coef1_t,coef2_t)));
-        const double lambda_l = std::max(0., std::min(1., std::min(coef1_l,coef2_l)));
-
-        auto coefficients = DPk_fh.coefficients(cell_id);
-        coefficients[0] = (1 - lambda_n)*fn*n[cell_id] + lambda_n*dot(coefficients[0],n[cell_id])*n[cell_id] 
-                        + (1 - lambda_t)*ft*t[cell_id] + lambda_t*dot(coefficients[0],t[cell_id])*t[cell_id] 
-                        + (1 - lambda_l)*fl*l[cell_id] + lambda_l*dot(coefficients[0],l[cell_id])*l[cell_id];
-        for(size_t i = 1; i < coefficients.size(); ++i){
-          coefficients[i] = lambda_n*dot(coefficients[i],n[cell_id])*n[cell_id] 
-                          + lambda_t*dot(coefficients[i],t[cell_id])*t[cell_id] 
-                          + lambda_l*dot(coefficients[i],l[cell_id])*l[cell_id];
-        }
-    });
-  }
-
-  void
-  _matrix_limiter(const MeshType& mesh,
-                  const DiscreteFunctionP0<const Rdxd>& S,
-                  DiscreteFunctionDPk<Dimension, Rdxd>& DPk_Sh) const
-  {
-    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
-    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
-    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
-    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
-                                                                        symmetry_boundary_descriptor_list);
-
-    const auto xj = mesh_data.xj();
-
-    //A retirer + tard
-    //const auto xr = mesh.xr();
-    //const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-    //
-
-    CellValue<double> lambda{mesh.connectivity()};
-    const DiscreteScalarFunction& inv2 = 0.5*trace(transpose(S)*S);
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-        const double inv2j = inv2[cell_id];
-        
-        double inv2_min = inv2j;
-        double inv2_max = inv2j;
-
-        const auto cell_stencil = stencil[cell_id];
-        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-          inv2_min = std::min(inv2_min, inv2[cell_stencil[i_cell]]);
-          inv2_max = std::max(inv2_max, inv2[cell_stencil[i_cell]]);
-        }
-
-        double inv2_bar_min = inv2j;
-        double inv2_bar_max = inv2j;
-
-        //const auto& cell_nodes = cell_to_node_matrix[cell_id];
-        //for(size_t r = 0; r < cell_nodes.size(); ++r){
-        //  const double inv2_xr = 0.5*trace(transpose(DPk_Sh[cell_id](xr[cell_nodes[r]]))*DPk_Sh[cell_id](xr[cell_nodes[r]]));
-//
-        //  inv2_bar_min = std::min(inv2_bar_min, inv2_xr);
-        //  inv2_bar_max = std::max(inv2_bar_max, inv2_xr);
-        //}
-
-        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-          const CellId cell_k_id = cell_stencil[i_cell];
-          const double inv2_xk   = 0.5*trace(transpose(DPk_Sh[cell_id](xj[cell_k_id]))*DPk_Sh[cell_id](xj[cell_k_id]));
-
-          inv2_bar_min = std::min(inv2_bar_min, inv2_xk);
-          inv2_bar_max = std::max(inv2_bar_max, inv2_xk);
-        }
-
-        const double eps =  1E-14;
-        double coef1     = 1.;
-        double coef2     = 1.;
-        if(std::abs(inv2_bar_max - inv2j) > eps){
-          coef1 = (inv2_max - inv2j) / (inv2_bar_max - inv2j);
-        } 
-        if(std::abs(inv2_bar_min - inv2j) > eps){
-          coef2 = (inv2_min - inv2j) / (inv2_bar_min - inv2j);
-        }
-
-        const double lambda_inv2 = sqrt(std::max(0., std::min(1., std::min(coef1, coef2))));
-        //const double lambda_inv2 = std::max(0., std::min(1., std::min(coef1, coef2)));
-        lambda[cell_id] = lambda_inv2;
-
-        auto coefficients = DPk_Sh.coefficients(cell_id);
-
-        coefficients[0] = (1-lambda_inv2)*S[cell_id] + lambda_inv2 * coefficients[0];
-        for (size_t i = 1; i < coefficients.size(); ++i) {
-          coefficients[i] *= lambda_inv2;
-        } 
-    });
-  } 
-
-  std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>>
-  _apply_NeoHook(const DiscreteScalarFunction rho,
-                 const DiscreteVectorFunction u,
-                 const DiscreteScalarFunction E,
-                 const DiscreteTensorFunction CG,
-                 const DiscreteScalarFunction lambda,
-                 const DiscreteScalarFunction mu,
-                 const DiscreteScalarFunction gamma,
-                 const DiscreteScalarFunction p_inf) const
-  {
-    const TinyMatrix<Dimension> I = identity;
-
-    const DiscreteScalarFunction& eps = E - 0.5 * dot(u, u);
-    const DiscreteScalarFunction& p_d = (gamma - 1) * rho * eps - p_inf * gamma;
-    const DiscreteTensorFunction& sigma_d =
-      (mu / sqrt(det(CG)) * (CG - I) + lambda / sqrt(det(CG)) * log(sqrt(det(CG))) * I) - p_d * I;
-    const DiscreteScalarFunction& aL_d = sqrt((lambda + 2 * mu) / rho + gamma * (p_d + p_inf) / rho);
-    const DiscreteScalarFunction& aT_d = sqrt(mu / rho);
-
-    return {std::make_shared<DiscreteFunctionVariant>(sigma_d), 
-            std::make_shared<DiscreteFunctionVariant>(p_d),
-            std::make_shared<DiscreteFunctionVariant>(aL_d),
-            std::make_shared<DiscreteFunctionVariant>(aT_d)};
-  }
-
-  std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>>
-  _apply_NeoHook(const DiscreteScalarFunction rho,
-                 const DiscreteVectorFunction u,
-                 const DiscreteScalarFunction E,
-                 const DiscreteTensorFunction CG,
-                 const DiscreteScalarFunction lambda,
-                 const DiscreteScalarFunction mu) const
-  {
-    const TinyMatrix<Dimension> I = identity;
-
-    const DiscreteScalarFunction& eps = E - 0.5 * dot(u, u);
-    const DiscreteTensorFunction& sigma_d =
-      (mu / sqrt(det(CG)) * (CG - I) + lambda / sqrt(det(CG)) * log(sqrt(det(CG))) * I);
-    const DiscreteScalarFunction& p_d = -1./3. * trace(sigma_d);
-    const DiscreteScalarFunction& aL_d = sqrt((lambda + 2 * mu) / rho);
-    const DiscreteScalarFunction& aT_d = sqrt(mu / rho);
-
-    return {std::make_shared<DiscreteFunctionVariant>(sigma_d), 
-            std::make_shared<DiscreteFunctionVariant>(p_d),
-            std::make_shared<DiscreteFunctionVariant>(aL_d),
-            std::make_shared<DiscreteFunctionVariant>(aT_d)};
-  }
-
-  std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>>
-  _apply_NeoHook2(const DiscreteScalarFunction rho,
-                  const DiscreteVectorFunction u,
-                  const DiscreteScalarFunction E,
-                  const DiscreteTensorFunction CG,
-                  const DiscreteScalarFunction mu,
-                  const DiscreteScalarFunction gamma,
-                  const DiscreteScalarFunction p_inf) const
-  {
-    const TinyMatrix<Dimension> I = identity;
-
-    const DiscreteScalarFunction& eps = E - 0.5 * dot(u, u);
-    const DiscreteScalarFunction& p_d =
-      (gamma - 1) * rho * eps - gamma * p_inf;
-    const DiscreteTensorFunction& sigma_d =
-      mu / sqrt(det(CG)) * (CG - 1. / 3. * trace(CG) * I) - p_d * I;
-    const DiscreteScalarFunction& aL_d = sqrt(
-      (2 * mu) / rho + (gamma * (p_d + p_inf)) / rho);
-    const DiscreteScalarFunction& aT_d = sqrt(mu / rho);
-
-    return {std::make_shared<DiscreteFunctionVariant>(sigma_d), 
-            std::make_shared<DiscreteFunctionVariant>(p_d),
-            std::make_shared<DiscreteFunctionVariant>(aL_d),
-            std::make_shared<DiscreteFunctionVariant>(aT_d)};
-  }
-
-  std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>,
-             const std::shared_ptr<const DiscreteFunctionVariant>>
-  _apply_state_law(const size_t law,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& CG_v,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& lambda_v,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& mu_v,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& gamma_v,
-                   const std::shared_ptr<const DiscreteFunctionVariant>& p_inf_v) const
-  {
-    const DiscreteScalarFunction& rho_d       = rho_v->get<DiscreteScalarFunction>();
-    const DiscreteVectorFunction& u_d         = u_v->get<DiscreteVectorFunction>();
-    const DiscreteScalarFunction& E_d         = E_v->get<DiscreteScalarFunction>();
-    const DiscreteTensorFunction& CG_d        = CG_v->get<DiscreteTensorFunction>();
-    const DiscreteScalarFunction& lambda_d    = lambda_v->get<DiscreteScalarFunction>();
-    const DiscreteScalarFunction& mu_d        = mu_v->get<DiscreteScalarFunction>();
-    const DiscreteScalarFunction& gamma_d     = gamma_v->get<DiscreteScalarFunction>();
-    const DiscreteScalarFunction& p_inf_d     = p_inf_v->get<DiscreteScalarFunction>();
-
-    if (law == 1) {
-      return _apply_NeoHook(rho_d, u_d, E_d, CG_d, lambda_d, mu_d, gamma_d, p_inf_d);
-    } else if (law == 2){
-      return _apply_NeoHook2(rho_d, u_d, E_d, CG_d, mu_d, gamma_d, p_inf_d);
-    } else {
-      return _apply_NeoHook(rho_d, u_d, E_d, CG_d,lambda_d, mu_d);
-    }
-  }
-
  public:
   std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
   compute_fluxes(const SolverType& solver_type,
@@ -1044,9 +619,9 @@ void
     DiscreteFunctionDPk<Dimension, double> p_lim   = copy(DPk_ph);
 
     const CellValue<const Rdxd> grad_u = this->_computeGradU(solver_type,rho_v,aL_v,aT_v,u_v,sigma_v,bc_descriptor_list);
-    this->_vector_limiter(mesh,u,u_lim,grad_u);
-    this->_matrix_limiter(mesh, S, S_lim);
-    this->_scalar_limiter(mesh, p, p_lim);
+    vector_limiter(mesh,u,u_lim,grad_u);
+    matrix_limiter(mesh, S, S_lim);
+    scalar_limiter(mesh, p, p_lim);
 
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
 
@@ -1473,7 +1048,7 @@ void
     std::shared_ptr<const DiscreteFunctionVariant> p_inf_app = shallowCopy(m_app, p_inf);
 
     auto [sigma_app, p_app, aL_app, aT_app] =
-    _apply_state_law(law,rho_app,u_app,E_app,CG_app,lambda_app,mu_app,gamma_app,p_inf_app);
+    apply_state_law<MeshType>(law,rho_app,u_app,E_app,CG_app,lambda_app,mu_app,gamma_app,p_inf_app);
 
     auto [ur_app, Fjr_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app, p_app, bc_descriptor_list);
     auto [ur_o1_app, Fjr_o1_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app,bc_descriptor_list);
diff --git a/src/scheme/Order2Limiters.hpp b/src/scheme/Order2Limiters.hpp
new file mode 100644
index 000000000..49ef8dcdf
--- /dev/null
+++ b/src/scheme/Order2Limiters.hpp
@@ -0,0 +1,330 @@
+#ifndef ORDER2_LIMITERS_HPP
+#define ORDER2_LIMITERS_HPP
+
+#include <algebra/EigenvalueSolver.hpp>
+
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/StencilArray.hpp>
+#include <mesh/StencilManager.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <scheme/PolynomialReconstructionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Socket.hpp>
+
+template <MeshConcept MeshType>
+void
+scalar_limiter(const MeshType& mesh,
+               const DiscreteFunctionP0<const double>& f,
+               DiscreteFunctionDPk<MeshType::Dimension, double>& DPk_fh)
+{
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
+
+    const auto xj = mesh_data.xj();
+
+    CellValue<double> alph_J2{mesh.connectivity()};
+
+    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+      const double fj = f[cell_id];
+
+      double f_min = fj;
+      double f_max = fj;
+
+      const auto cell_stencil = stencil[cell_id];
+      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+        f_min = std::min(f_min, f[cell_stencil[i_cell]]);
+        f_max = std::max(f_max, f[cell_stencil[i_cell]]);
+      }
+
+      double f_bar_min = fj;
+      double f_bar_max = fj;
+
+      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+        const CellId cell_k_id = cell_stencil[i_cell];
+        const double f_xk    = DPk_fh[cell_id](xj[cell_k_id]);
+
+        f_bar_min = std::min(f_bar_min, f_xk);
+        f_bar_max = std::max(f_bar_max, f_xk);
+      }
+
+      const double eps = 1E-14;
+      double coef1     = 1;
+      if (std::abs(f_bar_max - fj) > eps) {
+        coef1 = (f_max - fj) / ((f_bar_max - fj));
+      }
+
+      double coef2 = 1.;
+      if (std::abs(f_bar_min - fj) > eps) {
+        coef2 = (f_min - fj) / ((f_bar_min - fj));
+      }
+
+      const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
+      alph_J2[cell_id] = lambda;
+
+      auto coefficients = DPk_fh.coefficients(cell_id);
+
+      coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
+      for (size_t i = 1; i < coefficients.size(); ++i) {
+        coefficients[i] *= lambda;
+      }
+    });
+  }
+
+  template <MeshConcept MeshType>
+  void
+  vector_limiter(const MeshType& mesh,
+                 const DiscreteFunctionP0<const TinyVector<MeshType::Dimension>>& f,
+                 DiscreteFunctionDPk<MeshType::Dimension, TinyVector<MeshType::Dimension>> DPk_fh,
+                 const CellValue<const TinyMatrix<MeshType::Dimension>>& grad_u) 
+  {
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
+                                                                        symmetry_boundary_descriptor_list);
+    
+    const auto xj = mesh_data.xj();
+
+    CellValue<TinyMatrix<MeshType::Dimension>> sym_grad_u{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId j){
+        sym_grad_u[j] = (grad_u[j] + transpose(grad_u[j]));
+    });
+
+    CellValue<SmallArray<double>> eigunvalues{mesh.connectivity()};
+    CellValue<std::vector<SmallVector<double>>> eigunvectors{mesh.connectivity()};
+    for(CellId j = 0; j < mesh.numberOfCells(); ++j){
+      EigenvalueSolver{}.computeForSymmetricMatrix(sym_grad_u[j],eigunvalues[j],eigunvectors[j]);
+    }
+
+    CellValue<TinyVector<MeshType::Dimension>> n{mesh.connectivity()};
+    CellValue<TinyVector<MeshType::Dimension>> t{mesh.connectivity()};
+    CellValue<TinyVector<MeshType::Dimension>> l{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        n[cell_id] = zero;
+        t[cell_id] = zero;
+        l[cell_id] = zero;
+        //if(frobeniusNorm(sym_grad_u[cell_id]) > 1E-5){
+          if constexpr(MeshType::Dimension == 2){
+            TinyVector<MeshType::Dimension> nj;
+            TinyVector<MeshType::Dimension> tj;
+            for(size_t i = 0; i < MeshType::Dimension; ++i){
+              nj[i] = eigunvectors[cell_id][0][i];
+              tj[i] = eigunvectors[cell_id][1][i];
+            }
+            n[cell_id] = nj;
+            t[cell_id] = tj;
+            if(l2Norm(nj) != 0){
+              n[cell_id] = (1/l2Norm(nj))*n[cell_id];
+            }
+            if(l2Norm(tj) != 0){
+              t[cell_id] = (1/l2Norm(tj))*t[cell_id];
+            }
+          } else {
+            if constexpr(MeshType::Dimension == 3){
+              TinyVector<MeshType::Dimension> nj;
+              TinyVector<MeshType::Dimension> tj;
+              TinyVector<MeshType::Dimension> lj;
+              for(size_t i = 0; i < MeshType::Dimension; ++i){
+                nj[i] = eigunvectors[cell_id][0][i];
+                tj[i] = eigunvectors[cell_id][1][i];
+                lj[i] = eigunvectors[cell_id][2][i];
+              }
+              n[cell_id] = nj;
+              t[cell_id] = tj;
+              l[cell_id] = lj;
+              if(l2Norm(nj) != 0){
+              n[cell_id] = (1/l2Norm(nj))*n[cell_id];
+              }
+              if(l2Norm(tj) != 0){
+                t[cell_id] = (1/l2Norm(tj))*t[cell_id];
+              }
+              if(l2Norm(lj) != 0){
+                l[cell_id] = (1/l2Norm(lj))*l[cell_id];
+              }
+            } else {
+                TinyVector<MeshType::Dimension> nj;
+                nj[0] = 1.;
+                n[cell_id] = nj;
+            }
+          }
+        //}
+    });
+    
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+        const double fn = dot(f[cell_id],n[cell_id]);
+        const double ft = dot(f[cell_id],t[cell_id]);
+        const double fl = dot(f[cell_id],l[cell_id]);
+
+        double fn_min = fn;
+        double fn_max = fn;
+        double ft_min = ft;
+        double ft_max = ft;
+        double fl_min = fl;
+        double fl_max = fl;
+
+        const auto cell_stencil = stencil[cell_id];
+        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+          const double fn_k = dot(f[cell_stencil[i_cell]],n[cell_id]);
+          fn_min            = std::min(fn_min,fn_k);
+          fn_max            = std::max(fn_max,fn_k);
+
+          const double ft_k = dot(f[cell_stencil[i_cell]],t[cell_id]);
+          ft_min            = std::min(ft_min,ft_k);
+          ft_max            = std::max(ft_max,ft_k);
+
+          const double fl_k = dot(f[cell_stencil[i_cell]],l[cell_id]);
+          fl_min            = std::min(fl_min,fl_k);
+          fl_max            = std::max(fl_max,fl_k);
+        }
+
+        double fn_bar_min = fn;
+        double fn_bar_max = fn;
+        double ft_bar_min = ft;
+        double ft_bar_max = ft;
+        double fl_bar_min = fl;
+        double fl_bar_max = fl;
+
+        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+          const CellId cell_k_id = cell_stencil[i_cell];
+          const TinyVector<MeshType::Dimension> f_xk          = DPk_fh[cell_id](xj[cell_k_id]);
+
+          const double fn_xk = dot(f_xk,n[cell_id]);
+          fn_bar_min         = std::min(fn_bar_min,fn_xk);
+          fn_bar_max         = std::max(fn_bar_max,fn_xk);
+
+          const double ft_xk = dot(f_xk,t[cell_id]);
+          ft_bar_min         = std::min(ft_bar_min,ft_xk);
+          ft_bar_max         = std::max(ft_bar_max,ft_xk);
+
+          const double fl_xk = dot(f_xk,l[cell_id]);
+          fl_bar_min         = std::min(fl_bar_min,fl_xk);
+          fl_bar_max         = std::max(fl_bar_max,fl_xk);
+        }
+
+        const double eps = 1E-14;
+        double coef1_n   = 1;
+        if(std::abs(fn_bar_max - fn) > eps){
+          coef1_n = (fn_max - fn) / (fn_bar_max - fn);
+        } 
+        double coef2_n = 1;
+        if(std::abs(fn_bar_min - fn) > eps){
+          coef2_n = (fn_min - fn) / (fn_bar_min - fn);
+        }
+        
+        double coef1_t = 1;
+        if(std::abs(ft_bar_max - ft) > eps){
+          coef1_t = (ft_max - ft) / (ft_bar_max - ft);
+        }
+        double coef2_t = 1;
+        if(std::abs(ft_bar_min - ft) > eps){
+          coef2_t = (ft_min - ft) / (ft_bar_min - ft);
+        }
+
+        double coef1_l = 1;
+        if(std::abs(fl_bar_max - fl) > eps){
+          coef1_l = (fl_max - fl) / (fl_bar_max - fl);
+        }
+        double coef2_l = 1;
+        if(std::abs(fl_bar_min - fl) > eps){
+          coef2_l = (fl_min - fl) / (fl_bar_min - fl);
+        }
+
+        const double lambda_n = std::max(0., std::min(1., std::min(coef1_n,coef2_n)));
+        const double lambda_t = std::max(0., std::min(1., std::min(coef1_t,coef2_t)));
+        const double lambda_l = std::max(0., std::min(1., std::min(coef1_l,coef2_l)));
+
+        auto coefficients = DPk_fh.coefficients(cell_id);
+        coefficients[0] = (1 - lambda_n)*fn*n[cell_id] + lambda_n*dot(coefficients[0],n[cell_id])*n[cell_id] 
+                        + (1 - lambda_t)*ft*t[cell_id] + lambda_t*dot(coefficients[0],t[cell_id])*t[cell_id] 
+                        + (1 - lambda_l)*fl*l[cell_id] + lambda_l*dot(coefficients[0],l[cell_id])*l[cell_id];
+        for(size_t i = 1; i < coefficients.size(); ++i){
+          coefficients[i] = lambda_n*dot(coefficients[i],n[cell_id])*n[cell_id] 
+                          + lambda_t*dot(coefficients[i],t[cell_id])*t[cell_id] 
+                          + lambda_l*dot(coefficients[i],l[cell_id])*l[cell_id];
+        }
+    });
+  }
+
+  template <MeshConcept MeshType>
+  void
+  matrix_limiter(const MeshType& mesh,
+                 const DiscreteFunctionP0<const TinyMatrix<MeshType::Dimension>>& S,
+                 DiscreteFunctionDPk<MeshType::Dimension, TinyMatrix<MeshType::Dimension>>& DPk_Sh)
+  {
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
+                                                                        symmetry_boundary_descriptor_list);
+
+    const auto xj = mesh_data.xj();
+
+    const DiscreteFunctionP0<const double>& inv2 = 0.5*trace(transpose(S)*S);
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+        const double inv2j = inv2[cell_id];
+        
+        double inv2_min = inv2j;
+        double inv2_max = inv2j;
+
+        const auto cell_stencil = stencil[cell_id];
+        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+          inv2_min = std::min(inv2_min, inv2[cell_stencil[i_cell]]);
+          inv2_max = std::max(inv2_max, inv2[cell_stencil[i_cell]]);
+        }
+
+        double inv2_bar_min = inv2j;
+        double inv2_bar_max = inv2j;
+
+        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+          const CellId cell_k_id = cell_stencil[i_cell];
+          const double inv2_xk   = 0.5*trace(transpose(DPk_Sh[cell_id](xj[cell_k_id]))*DPk_Sh[cell_id](xj[cell_k_id]));
+
+          inv2_bar_min = std::min(inv2_bar_min, inv2_xk);
+          inv2_bar_max = std::max(inv2_bar_max, inv2_xk);
+        }
+
+        const double eps =  1E-14;
+        double coef1     = 1.;
+        double coef2     = 1.;
+        if(std::abs(inv2_bar_max - inv2j) > eps){
+          coef1 = (inv2_max - inv2j) / (inv2_bar_max - inv2j);
+        } 
+        if(std::abs(inv2_bar_min - inv2j) > eps){
+          coef2 = (inv2_min - inv2j) / (inv2_bar_min - inv2j);
+        }
+
+        const double lambda_inv2 = sqrt(std::max(0., std::min(1., std::min(coef1, coef2))));
+
+        auto coefficients = DPk_Sh.coefficients(cell_id);
+
+        coefficients[0] = (1-lambda_inv2)*S[cell_id] + lambda_inv2 * coefficients[0];
+        for (size_t i = 1; i < coefficients.size(); ++i) {
+          coefficients[i] *= lambda_inv2;
+        } 
+    });
+  } 
+
+#endif //ORDER2_LIMITERS_HPP
\ No newline at end of file
-- 
GitLab