#include <algebra/EigenvalueSolver.hpp>
#include <scheme/Order2LocalDtHyperelasticSolver.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/MeshVariant.hpp>
#include <mesh/StencilArray.hpp>
#include <mesh/StencilManager.hpp>
#include <mesh/SubItemValuePerItemVariant.hpp>
#include <scheme/CouplingBoundaryConditionDescriptor.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/HyperelasticSolver.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>
#include <utils/Socket.hpp>

#include <variant>
#include <vector>

template <MeshConcept MeshType>
class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver final
  : public Order2LocalDtHyperelasticSolverHandler::IOrder2LocalDtHyperelasticSolver
{
 private:
  static constexpr size_t Dimension = MeshType::Dimension;

  using Rdxd = TinyMatrix<Dimension>;
  using Rd   = TinyVector<Dimension>;

  using MeshDataType = MeshData<MeshType>;

  using DiscreteScalarFunction = DiscreteFunctionP0<const double>;
  using DiscreteVectorFunction = DiscreteFunctionP0<const Rd>;
  using DiscreteTensorFunction = DiscreteFunctionP0<const Rdxd>;

  class FixedBoundaryCondition;
  class PressureBoundaryCondition;
  class NormalStressBoundaryCondition;
  class SymmetryBoundaryCondition;
  class VelocityBoundaryCondition;
  class CouplingBoundaryCondition;

  using BoundaryCondition = std::variant<FixedBoundaryCondition,
                                         PressureBoundaryCondition,
                                         NormalStressBoundaryCondition,
                                         SymmetryBoundaryCondition,
                                         VelocityBoundaryCondition,
                                         CouplingBoundaryCondition>;

  using BoundaryConditionList = std::vector<BoundaryCondition>;

  NodeValuePerCell<const Rdxd>
  _computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoaL, const DiscreteScalarFunction& rhoaT) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
    const NodeValuePerCell<const Rd> njr = mesh_data.njr();

    NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
    const Rdxd I = identity;
    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
        const size_t& nb_nodes = Ajr.numberOfSubValues(j);
        const double& rhoaL_j  = rhoaL[j];
        const double& rhoaT_j  = rhoaT[j];
        for (size_t r = 0; r < nb_nodes; ++r) {
          const Rdxd& M = tensorProduct(Cjr(j, r), njr(j, r));
          Ajr(j, r)     = rhoaL_j * M + rhoaT_j * (l2Norm(Cjr(j, r)) * I - M);
        }
      });

    return Ajr;
  }

  NodeValuePerCell<const Rdxd>
  _computeEucclhydAjr(const MeshType& mesh,
                      const DiscreteScalarFunction& rhoaL,
                      const DiscreteScalarFunction& rhoaT) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
    const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();

    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
    const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();

    NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};

    parallel_for(
      Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Ajr[jr] = zero; });
    const Rdxd I = identity;
    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
        const auto& cell_nodes = cell_to_node_matrix[j];

        const auto& cell_faces = cell_to_face_matrix[j];

        const double& rho_aL = rhoaL[j];
        const double& rho_aT = rhoaT[j];

        for (size_t L = 0; L < cell_faces.size(); ++L) {
          const FaceId& l        = cell_faces[L];
          const auto& face_nodes = face_to_node_matrix[l];

          auto local_node_number_in_cell = [&](NodeId node_number) {
            for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
              if (node_number == cell_nodes[i_node]) {
                return i_node;
              }
            }
            return std::numeric_limits<size_t>::max();
          };

          for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
            const size_t R = local_node_number_in_cell(face_nodes[rl]);
            const Rdxd& M  = tensorProduct(Nlr(l, rl), nlr(l, rl));
            Ajr(j, R) += rho_aL * M + rho_aT * (l2Norm(Nlr(l, rl)) * I - M);
          }
        }
      });

    return Ajr;
  }

  NodeValuePerCell<const Rdxd>
  _computeAjr(const SolverType& solver_type,
              const MeshType& mesh,
              const DiscreteScalarFunction& rhoaL,
              const DiscreteScalarFunction& rhoaT) const
  {
    if constexpr (Dimension == 1) {
      return _computeGlaceAjr(mesh, rhoaL, rhoaT);
    } else {
      switch (solver_type) {
      case SolverType::Glace: {
        return _computeGlaceAjr(mesh, rhoaL, rhoaT);
      }
      case SolverType::Eucclhyd: {
        return _computeEucclhydAjr(mesh, rhoaL, rhoaT);
      }
      default: {
        throw UnexpectedError("invalid solver type");
      }
      }
    }
  }

  NodeValue<Rdxd>
  _computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr) const
  {
    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

    NodeValue<Rdxd> Ar{mesh.connectivity()};

    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
        Rdxd sum                                   = zero;
        const auto& node_to_cell                   = node_to_cell_matrix[r];
        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);

        for (size_t j = 0; j < node_to_cell.size(); ++j) {
          const CellId J       = node_to_cell[j];
          const unsigned int R = node_local_number_in_its_cells[j];
          sum += Ajr(J, R);
        }
        Ar[r] = sum;
      });

    return Ar;
  }

  NodeValue<Rd>
  _computeBr(const Mesh<Dimension>& mesh,
             const NodeValuePerCell<const Rdxd>& Ajr,
             DiscreteFunctionDPk<Dimension, const Rd> DPk_u,
             NodeValuePerCell<const Rdxd> DPk_sigma) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
    const auto& xr = mesh.xr();

    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

    NodeValue<Rd> b{mesh.connectivity()};

    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
        const auto& node_to_cell                   = node_to_cell_matrix[r];
        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);

        Rd br = zero;
        for (size_t j = 0; j < node_to_cell.size(); ++j) {
          const CellId J       = node_to_cell[j];
          const unsigned int R = node_local_number_in_its_cells[j];
          br += Ajr(J, R) * DPk_u[J](xr[r]) - DPk_sigma(J, R) * Cjr(J, R);
        }

        b[r] = br;
      });

    return b;
  }

  NodeValue<Rd>
  _computeBr(const Mesh<Dimension>& mesh,
             const NodeValuePerCell<const Rdxd>& Ajr,
             DiscreteFunctionDPk<Dimension, const Rd> DPk_u,
             DiscreteFunctionDPk<Dimension, const Rdxd> DPk_S,
             DiscreteFunctionDPk<Dimension, const double> DPk_p) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
    const auto& xr                        = mesh.xr();
    const TinyMatrix<Dimension> I = identity;

    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

    NodeValue<Rd> b{mesh.connectivity()};

    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
        const auto& node_to_cell                   = node_to_cell_matrix[r];
        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);

        Rd br = zero;
        for (size_t j = 0; j < node_to_cell.size(); ++j) {
          const CellId J       = node_to_cell[j];
          const unsigned int R = node_local_number_in_its_cells[j];
          const Rdxd sigma_r = DPk_S[J](xr[r]) - DPk_p[J](xr[r])*I;

          br += Ajr(J, R) * DPk_u[J](xr[r]) - sigma_r * Cjr(J, R);
        }

        b[r] = br;
      });

    return b;
  }

    NodeValue<Rd>
  _computeBr(const Mesh<Dimension>& mesh,
             const NodeValuePerCell<const Rdxd>& Ajr,
             const DiscreteVectorFunction& u,
             const DiscreteTensorFunction& sigma) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();

    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

    NodeValue<Rd> b{mesh.connectivity()};

    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
        const auto& node_to_cell                   = node_to_cell_matrix[r];
        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);

        Rd br = zero;
        for (size_t j = 0; j < node_to_cell.size(); ++j) {
          const CellId J       = node_to_cell[j];
          const unsigned int R = node_local_number_in_its_cells[j];
          br += Ajr(J, R) * u[J] - sigma[J] * Cjr(J, R);
        }

        b[r] = br;
      });

    return b;
  }

  BoundaryConditionList
  _getBCList(const MeshType& mesh,
             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
  {
    BoundaryConditionList bc_list;

    for (const auto& bc_descriptor : bc_descriptor_list) {
      bool is_valid_boundary_condition = true;

      switch (bc_descriptor->type()) {
      case IBoundaryConditionDescriptor::Type::symmetry: {
        bc_list.emplace_back(
          SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
        break;
      }
      case IBoundaryConditionDescriptor::Type::fixed: {
        bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
        break;
      }
      case IBoundaryConditionDescriptor::Type::dirichlet: {
        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
        if (dirichlet_bc_descriptor.name() == "velocity") {
          MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());

          Array<const Rd> value_list =
            InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
                                                                               mesh.xr(),
                                                                               mesh_node_boundary.nodeList());

          bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, value_list});
        } else if (dirichlet_bc_descriptor.name() == "pressure") {
          const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();

          if constexpr (Dimension == 1) {
            MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());

            Array<const double> node_values =
              InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(),
                                                                                     mesh_node_boundary.nodeList());

            bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
          } else {
            MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());

            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
            Array<const double> face_values =
              InterpolateItemValue<double(Rd)>::template interpolate<ItemType::face>(pressure_id, mesh_data.xl(),
                                                                                     mesh_face_boundary.faceList());
            bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values});
          }

        } else if (dirichlet_bc_descriptor.name() == "normal-stress") {
          const FunctionSymbolId normal_stress_id = dirichlet_bc_descriptor.rhsSymbolId();

          if constexpr (Dimension == 1) {
            MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());

            Array<const Rdxd> node_values =
              InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::node>(normal_stress_id, mesh.xr(),
                                                                                   mesh_node_boundary.nodeList());

            bc_list.emplace_back(NormalStressBoundaryCondition{mesh_node_boundary, node_values});
          } else {
            MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());

            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
            Array<const Rdxd> face_values =
              InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::face>(normal_stress_id, mesh_data.xl(),
                                                                                   mesh_face_boundary.faceList());
            bc_list.emplace_back(NormalStressBoundaryCondition{mesh_face_boundary, face_values});
          }

        } else {
          is_valid_boundary_condition = false;
        }
        break;
      }
      case IBoundaryConditionDescriptor::Type::coupling: {
        const CouplingBoundaryConditionDescriptor& coupling_bc_descriptor =
          dynamic_cast<const CouplingBoundaryConditionDescriptor&>(*bc_descriptor);
        bc_list.emplace_back(CouplingBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                                       coupling_bc_descriptor.label()));
        break;
      }
      default: {
        is_valid_boundary_condition = false;
      }
      }
      if (not is_valid_boundary_condition) {
        std::ostringstream error_msg;
        error_msg << *bc_descriptor << " is an invalid boundary condition for hyperelastic solver";
        throw NormalError(error_msg.str());
      }
    }

    return bc_list;
  }

  void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
  void _applyNormalStressBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
  void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
  void _applyVelocityBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
  void _applyCouplingBC(NodeValue<Rdxd>& Ar1,
                        NodeValue<Rdxd>& Ar2,
                        NodeValue<Rd>& br1,
                        NodeValue<Rd>& br2,
                        const std::vector<NodeId>& map1,
                        const std::vector<NodeId>& map2) const;
  void _applyCouplingBC(const MeshType& mesh,
                        NodeValue<Rd>& ur,
                        NodeValue<Rd>& CR_ur,
                        NodeValuePerCell<Rd>& Fjr,
                        NodeValuePerCell<Rd>& CR_Fjr,
                        const std::vector<NodeId>& map2) const;
  void _applyCouplingBC2(NodeValue<Rdxd>& Ar1,
                         NodeValue<Rdxd>& Ar2,
                         NodeValue<Rd>& ur1,
                         NodeValue<Rd>& ur2,
                         const std::vector<NodeId>& map1,
                         const std::vector<NodeId>& map2) const;
  void
  _applyBoundaryConditions(const BoundaryConditionList& bc_list,
                           const MeshType& mesh,
                           NodeValue<Rdxd>& Ar,
                           NodeValue<Rd>& br) const
  {
    this->_applyPressureBC(bc_list, mesh, br);
    this->_applyNormalStressBC(bc_list, mesh, br);
    this->_applySymmetryBC(bc_list, Ar, br);
    this->_applyVelocityBC(bc_list, Ar, br);
  }

  CellValue<Rdxd>
  _computeGradU(const SolverType& solver_type,
               const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
               const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
               const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
               const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
               const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
  {
    std::shared_ptr mesh_v = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
    if (not mesh_v) {
      throw NormalError("discrete functions are not defined on the same mesh");
    }

    if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
      throw NormalError("hyperelastic solver expects P0 functions");
    }

    const MeshType& mesh                = *mesh_v->get<MeshType>();
    const DiscreteScalarFunction& rho   = rho_v->get<DiscreteScalarFunction>();
    const DiscreteVectorFunction& u     = u_v->get<DiscreteVectorFunction>();
    const DiscreteScalarFunction& aL    = aL_v->get<DiscreteScalarFunction>();
    const DiscreteScalarFunction& aT    = aT_v->get<DiscreteScalarFunction>();
    const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();

    NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);

    NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, sigma);

    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
    this->_applyBoundaryConditions(bc_list, mesh, Ar, br);

    synchronize(Ar);
    synchronize(br);

    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();

    const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
    const NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);

    CellValue<Rdxd> grad_u{mesh.connectivity()};
    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
        auto cell_nodes = cell_to_node_matrix[j];

        Rdxd gradv           = zero;
        for (size_t R = 0; R < cell_nodes.size(); ++R) {
          NodeId r = cell_nodes[R];
          gradv += tensorProduct(ur[r], Cjr(j, R));
        }
        grad_u[j] = gradv;
      });

    return grad_u;
  }

  NodeValue<Rd>
  _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
  {
    NodeValue<Rd> u{mesh.connectivity()};
    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { u[r] = inverse(Ar[r]) * br[r]; });

    return u;
  }

  NodeValuePerCell<Rd>
  _computeFjr(const MeshType& mesh,
              const NodeValuePerCell<const Rdxd>& Ajr,
              const NodeValue<const Rd>& ur,
              DiscreteFunctionDPk<Dimension,const Rd> DPk_uh,
              NodeValuePerCell<const Rdxd> DPk_sigma) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
    const NodeValue<const Rd>& xr = mesh.xr();

    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

    NodeValuePerCell<Rd> F{mesh.connectivity()};
    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
        const auto& node_to_cell                   = node_to_cell_matrix[r];
        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);

        for(size_t j = 0; j < node_to_cell.size(); ++j){
          const CellId J       = node_to_cell[j];
          const unsigned int R = node_local_number_in_its_cells[j];
          F(J,R) = -Ajr(J,R) * (DPk_uh[J](xr[r]) - ur[r]) + DPk_sigma(J,R) * Cjr(J,R);
        }
      });

    return F;
  }

  NodeValuePerCell<Rd>
  _computeFjr(const MeshType& mesh,
              const NodeValuePerCell<const Rdxd>& Ajr,
              const NodeValue<const Rd>& ur,
              DiscreteFunctionDPk<Dimension, const Rd> DPk_uh,
              DiscreteFunctionDPk<Dimension, const Rdxd> DPk_S,
              DiscreteFunctionDPk<Dimension, const double> DPk_p) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
    const NodeValue<const Rd>& xr        = mesh.xr();
    const TinyMatrix<Dimension> I = identity;

    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();

    NodeValuePerCell<Rd> F{mesh.connectivity()};

    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
        const auto& cell_nodes = cell_to_node_matrix[j];

        for(size_t r = 0; r < cell_nodes.size(); ++r) {
          Rdxd sigma_r = DPk_S[j](xr[cell_nodes[r]]) - DPk_p[j](xr[cell_nodes[r]])*I;
          F(j,r) = -Ajr(j,r) * (DPk_uh[j](xr[cell_nodes[r]]) - ur[cell_nodes[r]]) + sigma_r*Cjr(j,r);
        }
      });

    return F;
  }

  NodeValuePerCell<Rd>
  _computeFjr(const MeshType& mesh,
              const NodeValuePerCell<const Rdxd>& Ajr,
              const NodeValue<const Rd>& ur,
              const DiscreteVectorFunction& u,
              const DiscreteTensorFunction& sigma) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();

    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();

    NodeValuePerCell<Rd> F{mesh.connectivity()};
    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
        const auto& cell_nodes = cell_to_node_matrix[j];

        for (size_t r = 0; r < cell_nodes.size(); ++r) {
          F(j, r) = -Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + sigma[j] * Cjr(j, r);
        }
      });

    return F;
  }

  std::tuple<std::vector<NodeId>, std::vector<NodeId>>
  _computeMapping(std::shared_ptr<const MeshVariant>& i_mesh1,
                  std::shared_ptr<const MeshVariant>& i_mesh2,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
  {
    const MeshType& mesh1 = *i_mesh1->get<MeshType>();
    const MeshType& mesh2 = *i_mesh2->get<MeshType>();

    const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
    const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);

    std::vector<NodeId> map1;
    std::vector<NodeId> map2;

    NodeValue<Rd> xr1 = copy(mesh1.xr());
    NodeValue<Rd> xr2 = copy(mesh2.xr());
    for (const auto& boundary_condition1 : bc_list1) {
      std::visit(
        [&](auto&& bc1) {
          using T1 = std::decay_t<decltype(bc1)>;
          if constexpr (std::is_same_v<CouplingBoundaryCondition, T1>) {
            const auto& node_list1 = bc1.nodeList();
            for (const auto& boundary_condition2 : bc_list2) {
              std::visit(
                [&](auto&& bc2) {
                  using T2 = std::decay_t<decltype(bc2)>;
                  if constexpr (std::is_same_v<CouplingBoundaryCondition, T2>) {
                    if (bc1.label() == bc2.label()) {
                      const auto& node_list2 = bc2.nodeList();
                      for (size_t i = 0; i < node_list1.size(); i++) {
                        NodeId node_id1 = node_list1[i];
                        NodeId node_id2;
                        for (size_t j = 0; j < node_list2.size(); j++) {
                          node_id2   = node_list2[j];
                          double err = 0;
                          err        = dot((xr1[node_id1] - xr2[node_id2]), (xr1[node_id1] - xr2[node_id2]));
                          if (sqrt(err) < 1e-14) {
                            map1.push_back(node_id1);
                            map2.push_back(node_id2);
                          }
                        }
                      }
                    }
                  }
                },
                boundary_condition2);
            }
          }
        },
        boundary_condition1);
    }
    return std::make_tuple(map1, map2);
  }

 public:
  std::tuple<const std::shared_ptr<const ItemValueVariant>,
             const std::shared_ptr<const SubItemValuePerItemVariant>,
             const std::shared_ptr<const ItemValueVariant>,
             const std::shared_ptr<const SubItemValuePerItemVariant>,
             NodeValue<Rd>,
             NodeValuePerCell<Rd>>
  compute_fluxes(const SolverType& solver_type,
                 const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
                 const std::vector<NodeId>& map1,
                 const std::vector<NodeId>& map2) const
  {
    std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, aL1_v, aT1_v, u1_v, sigma1_v});
    std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, aL2_v, aT2_v, u2_v, sigma2_v});
    if (not i_mesh1) {
      throw NormalError("discrete functions are not defined on the same mesh");
    }

    if (not checkDiscretizationType({rho1_v, aL1_v, u1_v, sigma1_v}, DiscreteFunctionType::P0)) {
      throw NormalError("hyperelastic solver expects P0 functions");
    }
    if (not i_mesh2) {
      throw NormalError("discrete functions are not defined on the same mesh");
    }

    if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
      throw NormalError("acoustic solver expects P0 functions");
    }

    const MeshType& mesh1                = *i_mesh1->get<MeshType>();
    const DiscreteScalarFunction& rho1   = rho1_v->get<DiscreteScalarFunction>();
    const DiscreteVectorFunction& u1     = u1_v->get<DiscreteVectorFunction>();
    const DiscreteScalarFunction& aL1    = aL1_v->get<DiscreteScalarFunction>();
    const DiscreteScalarFunction& aT1    = aT1_v->get<DiscreteScalarFunction>();
    const DiscreteTensorFunction& sigma1 = sigma1_v->get<DiscreteTensorFunction>();
    const DiscreteScalarFunction& p1     = p1_v->get<DiscreteScalarFunction>();

    const MeshType& mesh2                = *i_mesh2->get<MeshType>();
    const DiscreteScalarFunction& rho2   = rho2_v->get<DiscreteScalarFunction>();
    const DiscreteVectorFunction& u2     = u2_v->get<DiscreteVectorFunction>();
    const DiscreteScalarFunction& aL2    = aL2_v->get<DiscreteScalarFunction>();
    const DiscreteScalarFunction& aT2    = aT2_v->get<DiscreteScalarFunction>();
    const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
    const DiscreteScalarFunction& p2     = p2_v->get<DiscreteScalarFunction>();

    const TinyMatrix<Dimension> I    = identity;
    const DiscreteTensorFunction& S1 = sigma1 + p1*I;
    const DiscreteTensorFunction& S2 = sigma2 + p2*I;

    const std::shared_ptr<const DiscreteFunctionVariant>& S1_v = std::make_shared<DiscreteFunctionVariant>(S1);
    const std::shared_ptr<const DiscreteFunctionVariant>& S2_v = std::make_shared<DiscreteFunctionVariant>(S2);

    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list1;
    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list2;

    for(auto&& bc_descriptor : bc_descriptor_list1){
      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
        symmetry_boundary_descriptor_list1.push_back(bc_descriptor->boundaryDescriptor_shared());
      }
    }
    for(auto&& bc_descriptor : bc_descriptor_list2){
      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
        symmetry_boundary_descriptor_list2.push_back(bc_descriptor->boundaryDescriptor_shared());
      }
    }

    PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1,
                                                                 symmetry_boundary_descriptor_list1);
    PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1,
                                                                 symmetry_boundary_descriptor_list2);

    auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v, S1_v, p1_v);
    auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
    auto DPk_Sh1 = reconstruction1[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
    auto DPk_ph1 = reconstruction1[2]->get<DiscreteFunctionDPk<Dimension, const double>>();

    DiscreteFunctionDPk<Dimension, Rd> u1_lim       = copy(DPk_uh1);
    DiscreteFunctionDPk<Dimension, Rdxd> S1_lim     = copy(DPk_Sh1);
    DiscreteFunctionDPk<Dimension, double> p1_lim   = copy(DPk_ph1);

    const CellValue<const Rdxd> grad_u1 = this->_computeGradU(solver_type,rho1_v,aL1_v,aT1_v,u1_v,sigma1_v,bc_descriptor_list1);
    vector_limiter(mesh1,u1,u1_lim,grad_u1);
    matrix_limiter(mesh1, S1, S1_lim);
    scalar_limiter(mesh1, p1, p1_lim);

    auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v, S2_v, p2_v);
    auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
    auto DPk_Sh2 = reconstruction2[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
    auto DPk_ph2 = reconstruction2[2]->get<DiscreteFunctionDPk<Dimension, const double>>();

    DiscreteFunctionDPk<Dimension, Rd> u2_lim       = copy(DPk_uh2);
    DiscreteFunctionDPk<Dimension, Rdxd> S2_lim     = copy(DPk_Sh2);
    DiscreteFunctionDPk<Dimension, double> p2_lim   = copy(DPk_ph2);

    const CellValue<const Rdxd> grad_u2 = this->_computeGradU(solver_type,rho2_v,aL2_v,aT2_v,u2_v,sigma2_v,bc_descriptor_list2);
    vector_limiter(mesh2,u2,u2_lim,grad_u2);
    matrix_limiter(mesh2, S2, S2_lim);
    scalar_limiter(mesh2, p2, p2_lim);

    NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
    NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);

    NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1_lim, S1_lim, p1_lim);
    NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2_lim, S2_lim, p2_lim);

    const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
    const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
    this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
    this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);

    synchronize(Ar1);
    synchronize(br1);
    synchronize(Ar2);
    synchronize(br2);

    NodeValue<Rd> ur1           = this->_computeUr(mesh1, Ar1, br1);
    NodeValue<Rd> ur2           = this->_computeUr(mesh2, Ar2, br2);

    this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);

    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, S1_lim, p1_lim);
    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, S2_lim, p2_lim);

    NodeValue<Rd> CR_ur         = copy(ur2);
    NodeValuePerCell<Rd> CR_Fjr = copy(Fjr2);

    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
                           std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
                           std::make_shared<const ItemValueVariant>(ur2),
                           std::make_shared<const SubItemValuePerItemVariant>(Fjr2), CR_ur, CR_Fjr);
  }

  std::tuple<const std::shared_ptr<const ItemValueVariant>,
             const std::shared_ptr<const SubItemValuePerItemVariant>,
             const std::shared_ptr<const ItemValueVariant>,
             const std::shared_ptr<const SubItemValuePerItemVariant>,
             NodeValue<Rd>,
             NodeValuePerCell<Rd>>
  compute_fluxes(const SolverType& solver_type,
                 const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_v,
                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2_v,
                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
                 const std::vector<NodeId>& map1,
                 const std::vector<NodeId>& map2) const
  {
    std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, aL1_v, aT1_v, u1_v, sigma1_v});
    std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, aL2_v, aT2_v, u2_v, sigma2_v});
    if (not i_mesh1) {
      throw NormalError("discrete functions are not defined on the same mesh");
    }

    if (not checkDiscretizationType({rho1_v, aL1_v, u1_v, sigma1_v}, DiscreteFunctionType::P0)) {
      throw NormalError("hyperelastic solver expects P0 functions");
    }
    if (not i_mesh2) {
      throw NormalError("discrete functions are not defined on the same mesh");
    }

    if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
      throw NormalError("acoustic solver expects P0 functions");
    }

    const MeshType& mesh1                = *i_mesh1->get<MeshType>();
    const DiscreteScalarFunction& rho1   = rho1_v->get<DiscreteScalarFunction>();
    const DiscreteVectorFunction& u1     = u1_v->get<DiscreteVectorFunction>();
    const DiscreteScalarFunction& aL1    = aL1_v->get<DiscreteScalarFunction>();
    const DiscreteScalarFunction& aT1    = aT1_v->get<DiscreteScalarFunction>();
    const DiscreteTensorFunction& sigma1 = sigma1_v->get<DiscreteTensorFunction>();

    const MeshType& mesh2                = *i_mesh2->get<MeshType>();
    const DiscreteScalarFunction& rho2   = rho2_v->get<DiscreteScalarFunction>();
    const DiscreteVectorFunction& u2     = u2_v->get<DiscreteVectorFunction>();
    const DiscreteScalarFunction& aL2    = aL2_v->get<DiscreteScalarFunction>();
    const DiscreteScalarFunction& aT2    = aT2_v->get<DiscreteScalarFunction>();
    const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();

    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list1;
    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list2;

    for(auto&& bc_descriptor : bc_descriptor_list1){
      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
        symmetry_boundary_descriptor_list1.push_back(bc_descriptor->boundaryDescriptor_shared());
      }
    }
    for(auto&& bc_descriptor : bc_descriptor_list2){
      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
        symmetry_boundary_descriptor_list2.push_back(bc_descriptor->boundaryDescriptor_shared());
      }
    }

    NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
    NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);

    NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1, sigma1);
    NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2, sigma2);

    const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
    const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
    this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
    this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);

    synchronize(Ar1);
    synchronize(br1);
    synchronize(Ar2);
    synchronize(br2);

    NodeValue<Rd> ur1           = this->_computeUr(mesh1, Ar1, br1);
    NodeValue<Rd> ur2           = this->_computeUr(mesh2, Ar2, br2);

    this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);

    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1);
    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);

    NodeValue<Rd> CR_ur         = copy(ur2);
    NodeValuePerCell<Rd> CR_Fjr = copy(Fjr2);

    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
                           std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
                           std::make_shared<const ItemValueVariant>(ur2),
                           std::make_shared<const SubItemValuePerItemVariant>(Fjr2), CR_ur, CR_Fjr);
  }

  std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
  compute_fluxes(const SolverType& solver_type,
                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& p_v, 
                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                 NodeValue<Rd> CR_ur,
                 NodeValuePerCell<Rd> CR_Fjr,
                 const std::vector<NodeId> map2) const
  {
    std::shared_ptr i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
    if (not i_mesh) {
      throw NormalError("discrete functions are not defined on the same mesh ici1");
    }

    if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
      throw NormalError("hyperelastic solver expects P0 functions");
    }

    const MeshType& mesh                = *i_mesh->get<MeshType>();
    const DiscreteScalarFunction& rho   = rho_v->get<DiscreteScalarFunction>();
    const DiscreteVectorFunction& u     = u_v->get<DiscreteVectorFunction>();
    const DiscreteScalarFunction& aL    = aL_v->get<DiscreteScalarFunction>();
    const DiscreteScalarFunction& aT    = aT_v->get<DiscreteScalarFunction>();
    const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
    const DiscreteScalarFunction& p     = p_v->get<DiscreteScalarFunction>();

    const TinyMatrix<Dimension> I = identity;
    const DiscreteTensorFunction& S     = sigma + p*I;
    const std::shared_ptr<const DiscreteFunctionVariant>& S_v = std::make_shared<DiscreteFunctionVariant>(S);

    std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;

    for(auto&& bc_descriptor : bc_descriptor_list){
      if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
        symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
      }
    }

    PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
                                                                symmetry_boundary_descriptor_list);

    auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v, S_v, p_v);
    auto DPk_uh     = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
    auto DPk_Sh     = reconstruction[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
    auto DPk_ph     = reconstruction[2]->get<DiscreteFunctionDPk<Dimension, const double>>();

    DiscreteFunctionDPk<Dimension, Rd> u_lim       = copy(DPk_uh);
    DiscreteFunctionDPk<Dimension, Rdxd> S_lim     = copy(DPk_Sh);
    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);
    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);

    NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u_lim, S_lim, p_lim);

    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
    this->_applyBoundaryConditions(bc_list, mesh, Ar, br);

    synchronize(Ar);
    synchronize(br);

    NodeValue<Rd> ur         = this->_computeUr(mesh, Ar, br);
    NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, S_lim, p_lim);

    this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);

    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                           std::make_shared<const SubItemValuePerItemVariant>(Fjr));
  }

  std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
  compute_fluxes(const SolverType& solver_type,
                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                 NodeValue<Rd> CR_ur,
                 NodeValuePerCell<Rd> CR_Fjr,
                 const std::vector<NodeId> map2) const
  {
    std::shared_ptr i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
    if (not i_mesh) {
      throw NormalError("discrete functions are not defined on the same mesh");
    }

    if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
      throw NormalError("hyperelastic solver expects P0 functions");
    }

    const MeshType& mesh                = *i_mesh->get<MeshType>();
    const DiscreteScalarFunction& rho   = rho_v->get<DiscreteScalarFunction>();
    const DiscreteVectorFunction& u     = u_v->get<DiscreteVectorFunction>();
    const DiscreteScalarFunction& aL    = aL_v->get<DiscreteScalarFunction>();
    const DiscreteScalarFunction& aT    = aT_v->get<DiscreteScalarFunction>();
    const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();

    NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);

    NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, sigma);

    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
    this->_applyBoundaryConditions(bc_list, mesh, Ar, br);

    synchronize(Ar);
    synchronize(br);

    NodeValue<Rd> ur         = this->_computeUr(mesh, Ar, br);
    NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);

    this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);

    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                           std::make_shared<const SubItemValuePerItemVariant>(Fjr));
  }

  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>>
  apply_fluxes(const double& dt,
               const MeshType& mesh,
               const DiscreteScalarFunction& rho,
               const DiscreteVectorFunction& u,
               const DiscreteScalarFunction& E,
               const DiscreteTensorFunction& CG,
               const NodeValue<const Rd>& ur_o2,
               const NodeValue<const Rd>& ur_o1,
               const NodeValuePerCell<const Rd>& Fjr_o2,
               const NodeValuePerCell<const Rd>& Fjr_o1) const
  {
    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();

    NodeValue<Rd> ur         = copy(ur_o2);
    NodeValuePerCell<Rd> Fjr = copy(Fjr_o2);

    if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
        (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
      throw NormalError("fluxes are not defined on the same connectivity than the 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);

    CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
    const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();

    CellValue<double> new_rho = copy(rho.cellValues());
    CellValue<Rd> new_u       = copy(u.cellValues());
    CellValue<double> new_E   = copy(E.cellValues());
    CellValue<Rdxd> new_CG    = copy(CG.cellValues());

    CellValue<double> new_rho_stencil = copy(rho.cellValues());
    CellValue<Rd> new_u_stencil       = copy(u.cellValues());
    CellValue<double> new_E_stencil   = copy(E.cellValues());
    CellValue<Rdxd> new_CG_stencil    = copy(CG.cellValues());

    CellValue<int> cell_flag{mesh.connectivity()};

    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
        auto cell_nodes = cell_to_node_matrix[j];

        Rd momentum_fluxes   = zero;
        double energy_fluxes = 0;
        Rdxd gradv           = zero;
        for (size_t R = 0; R < cell_nodes.size(); ++R) {
          NodeId r = cell_nodes[R];
          gradv += tensorProduct(ur[r], Cjr(j, R));
          momentum_fluxes += Fjr(j, R);
          energy_fluxes += dot(Fjr(j, R), ur[r]);
        }
        Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
        double dt_over_Mj        = dt / (rho[j] * Vj[j]);
        double dt_over_Vj        = dt / Vj[j];
        new_u[j] += dt_over_Mj * momentum_fluxes;
        new_E[j] += dt_over_Mj * energy_fluxes;
        new_CG[j] += dt_over_Vj * cauchy_green_fluxes;
        new_CG[j] += transpose(new_CG[j]);
        new_CG[j] *= 0.5;

        double new_eps = new_E[j] - 0.5 * dot(new_u[j],new_u[j]);
        if(new_eps < 0.){
          cell_flag[j] = 1;
        }else{
          cell_flag[j] = 0;
        }
      });

    for(CellId j = 0; j < mesh.numberOfCells(); ++j){
      if(cell_flag[j] == 1){
        const auto& cell_nodes = cell_to_node_matrix[j];

        Rd momentum_fluxes   = zero;
        double energy_fluxes = 0;
        Rdxd gradv           = zero;
        std::vector<NodeId> NodeList;
        for (size_t R = 0; R < cell_nodes.size(); ++R) {
          NodeId r = cell_nodes[R];
          NodeList.push_back(r);
          Fjr(j,R) = Fjr_o1(j,R);
          ur[r] = ur_o1[r];
          gradv += tensorProduct(ur[r], Cjr(j, R));
          momentum_fluxes += Fjr(j, R);
          energy_fluxes += dot(Fjr(j, R), ur[r]);
        }
        Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
        double dt_over_Mj        = dt / (rho[j] * Vj[j]);
        double dt_over_Vj        = dt / Vj[j];
        new_u[j] = new_u_stencil[j] +  dt_over_Mj * momentum_fluxes;
        new_E[j] = new_E_stencil[j] + dt_over_Mj * energy_fluxes;
        new_CG[j] = new_CG_stencil[j] + dt_over_Vj * cauchy_green_fluxes;
        new_CG[j] += transpose(new_CG[j]);
        new_CG[j] *= 0.5;

        auto cell_stencil = stencil[j];
        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
          CellId J = cell_stencil[i_cell];
          auto i_cell_nodes = cell_to_node_matrix[J];
          momentum_fluxes = zero;
          energy_fluxes   = 0;
          gradv           = zero;
          for (size_t R = 0; R < i_cell_nodes.size(); ++R) {
            NodeId i_node = i_cell_nodes[R];
            for(size_t node_test = 0; node_test < NodeList.size(); ++node_test){
              if(NodeList[node_test] == i_node){
                Fjr(J,R) = Fjr_o1(J,R);
              }
            }
            gradv += tensorProduct(ur[i_node], Cjr(J, R));
            momentum_fluxes += Fjr(J, R);
            energy_fluxes += dot(Fjr(J, R), ur[i_node]);
          }
          cauchy_green_fluxes = gradv * CG[J] + CG[J] * transpose(gradv);
          dt_over_Mj          = dt / (rho[J] * Vj[J]);
          dt_over_Vj          = dt / Vj[J];
          new_u[J] = new_u_stencil[J] +  dt_over_Mj * momentum_fluxes;
          new_E[J] = new_E_stencil[J] + dt_over_Mj * energy_fluxes;
          new_CG[J] = new_CG_stencil[J] + dt_over_Vj * cauchy_green_fluxes;
          new_CG[J] += transpose(new_CG[J]);
          new_CG[J] *= 0.5;
        } 
      }
    }

    NodeValue<Rd> new_xr = copy(mesh.xr());
    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });

    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);

    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();

    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });

    return {std::make_shared<MeshVariant>(new_mesh),
            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
            std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_CG))};
  }

  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>>
  apply_fluxes(const double& dt,
               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,
               const std::shared_ptr<const ItemValueVariant>& ur,
               const std::shared_ptr<const ItemValueVariant>& ur_o1,
               const std::shared_ptr<const SubItemValuePerItemVariant> Fjr,
               const std::shared_ptr<const SubItemValuePerItemVariant> Fjr_o1) const
  {
    std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
    if (not mesh_v) {
      throw NormalError("discrete functions are not defined on the same mesh");
    }

    if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
      throw NormalError("hyperelastic solver expects P0 functions");
    }

    return this->apply_fluxes(dt,                                   //
                              *mesh_v->get<MeshType>(),             //
                              rho->get<DiscreteScalarFunction>(),   //
                              u->get<DiscreteVectorFunction>(),     //
                              E->get<DiscreteScalarFunction>(),     //
                              CG->get<DiscreteTensorFunction>(),    //
                              ur->get<NodeValue<const Rd>>(),
                              ur_o1->get<NodeValue<const Rd>>(),       //
                              Fjr->get<NodeValuePerCell<const Rd>>(),
                              Fjr_o1->get<NodeValuePerCell<const Rd>>()
                              );
  }

  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>>
  order2_apply_fluxes(const double& dt,
                      const MeshType& mesh,
                      const DiscreteScalarFunction& rho,
                      const DiscreteVectorFunction& u,
                      const DiscreteScalarFunction& E,
                      const DiscreteTensorFunction& CG,
                      const DiscreteScalarFunction& rho_app,
                      const DiscreteTensorFunction& CG_app,
                      const NodeValue<const Rd>& ur_o2,
                      const NodeValue<const Rd>& ur_o1,
                      const NodeValuePerCell<const Rd>& Fjr_o2,
                      const NodeValuePerCell<const Rd>& Fjr_o1) const
  {
    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();

    NodeValue<Rd> ur         = copy(ur_o2);
    NodeValuePerCell<Rd> Fjr = copy(Fjr_o2);

    if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
        (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
      throw NormalError("fluxes are not defined on the same connectivity than the 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);

    CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
    const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();

    CellValue<double> new_rho = copy(rho.cellValues());
    CellValue<Rd> new_u       = copy(u.cellValues());
    CellValue<double> new_E   = copy(E.cellValues());
    CellValue<Rdxd> new_CG    = copy(CG.cellValues());

    CellValue<double> new_rho_stencil = copy(rho.cellValues());
    CellValue<Rd> new_u_stencil       = copy(u.cellValues());
    CellValue<double> new_E_stencil   = copy(E.cellValues());
    CellValue<Rdxd> new_CG_stencil    = copy(CG.cellValues());

    CellValue<int> cell_flag{mesh.connectivity()};

    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
        const auto& cell_nodes = cell_to_node_matrix[j];

        Rd momentum_fluxes   = zero;
        double energy_fluxes = 0;
        Rdxd gradv           = zero;
        for (size_t R = 0; R < cell_nodes.size(); ++R) {
          const NodeId r = cell_nodes[R];
          gradv += tensorProduct(ur[r], Cjr(j, R));
          momentum_fluxes += Fjr(j, R);
          energy_fluxes += dot(Fjr(j, R), ur[r]);
        }
        const Rdxd cauchy_green_fluxes = rho_app[j] * (gradv * CG_app[j] + CG_app[j] * transpose(gradv));
        const double dt_over_Mj        = dt / (rho[j] * Vj[j]);
        new_u[j] += dt_over_Mj * momentum_fluxes;
        new_E[j] += dt_over_Mj * energy_fluxes;
        new_CG[j] += dt_over_Mj * cauchy_green_fluxes;
        new_CG[j] += transpose(new_CG[j]);
        new_CG[j] *= 0.5;

        const double& new_eps = new_E[j] - 0.5 * dot(new_u[j],new_u[j]);
        if(new_eps < 0.){
          cell_flag[j] = 1;
        }else{
          cell_flag[j] = 0;
        }
      });

    for(CellId j = 0; j < mesh.numberOfCells(); ++j){
      if(cell_flag[j] == 1){
        const auto& cell_nodes = cell_to_node_matrix[j];

        Rd momentum_fluxes   = zero;
        double energy_fluxes = 0;
        Rdxd gradv           = zero;
        std::vector<NodeId> NodeList;
        for (size_t R = 0; R < cell_nodes.size(); ++R) {
          NodeId r = cell_nodes[R];
          NodeList.push_back(r);
          Fjr(j,R) = Fjr_o1(j,R);
          ur[r] = ur_o1[r];
          gradv += tensorProduct(ur[r], Cjr(j, R));
          momentum_fluxes += Fjr(j, R);
          energy_fluxes += dot(Fjr(j, R), ur[r]);
        }
        Rdxd cauchy_green_fluxes = rho_app[j] * (gradv * CG_app[j] + CG_app[j] * transpose(gradv));
        double dt_over_Mj        = dt / (rho[j] * Vj[j]);
        new_u[j] = new_u_stencil[j] + dt_over_Mj * momentum_fluxes;
        new_E[j] = new_E_stencil[j] + dt_over_Mj * energy_fluxes;
        new_CG[j] = new_CG_stencil[j] + dt_over_Mj * cauchy_green_fluxes;
        new_CG[j] += transpose(new_CG[j]);
        new_CG[j] *= 0.5;

        auto cell_stencil = stencil[j];
        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
          CellId J = cell_stencil[i_cell];
          auto i_cell_nodes = cell_to_node_matrix[J];
          momentum_fluxes = zero;
          energy_fluxes   = 0;
          gradv           = zero;
          for (size_t R = 0; R < i_cell_nodes.size(); ++R) {
            NodeId i_node = i_cell_nodes[R];
            for(size_t node_test = 0; node_test < NodeList.size(); ++node_test){
              if(NodeList[node_test] == i_node){
                Fjr(J,R) = Fjr_o1(J,R);
              }
            }
            gradv += tensorProduct(ur[i_node], Cjr(J, R));
            momentum_fluxes += Fjr(J, R);
            energy_fluxes += dot(Fjr(J, R), ur[i_node]);
          }
          cauchy_green_fluxes = rho_app[J] * (gradv * CG_app[J] + CG_app[J] * transpose(gradv));
          dt_over_Mj        = dt / (rho[J] * Vj[J]);
          new_u[J] = new_u_stencil[J] + dt_over_Mj * momentum_fluxes;
          new_E[J] = new_E_stencil[J] + dt_over_Mj * energy_fluxes;
          new_CG[J] = new_CG_stencil[J] + dt_over_Mj * cauchy_green_fluxes;
          new_CG[J] += transpose(new_CG[J]);
          new_CG[J] *= 0.5;
        } 
      }
    }

    NodeValue<Rd> new_xr = copy(mesh.xr());
    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });

    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);

    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();

    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });

    return {std::make_shared<MeshVariant>(new_mesh),
            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
            std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_CG))};
  }

  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>>
  order2_apply_fluxes(const double& dt,
                      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,
                      const std::shared_ptr<const DiscreteFunctionVariant>& rho_app,
                      const std::shared_ptr<const DiscreteFunctionVariant>& CG_app,
                      const std::shared_ptr<const ItemValueVariant>& ur,
                      const std::shared_ptr<const ItemValueVariant>& ur_o1,
                      const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
                      const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_o1) const
  {
    std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
    if (not mesh_v) {
      throw NormalError("discrete functions are not defined on the same mesh");
    }

    if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
      throw NormalError("hyperelastic solver expects P0 functions");
    }

    return this->order2_apply_fluxes(dt,                                       //
                                     *mesh_v->get<MeshType>(),                 //
                                     rho->get<DiscreteScalarFunction>(),       //
                                     u->get<DiscreteVectorFunction>(),         //
                                     E->get<DiscreteScalarFunction>(),         //
                                     CG->get<DiscreteTensorFunction>(),        //
                                     rho_app->get<DiscreteScalarFunction>(),   //
                                     CG_app->get<DiscreteTensorFunction>(),    //
                                     ur->get<NodeValue<const Rd>>(),           //
                                     ur_o1->get<NodeValue<const Rd>>(),
                                     Fjr->get<NodeValuePerCell<const Rd>>(),
                                     Fjr_o1->get<NodeValuePerCell<const Rd>>());
  }

  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>>
  compute_sub_iterations(const SolverType& solver_type,
                         const size_t& law,
                         const double& CFL,
                         const double& dt,
                         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,
                         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                         NodeValue<Rd> CR_ur,
                         NodeValuePerCell<Rd> CR_Fjr,
                         std::vector<NodeId> map2,
                         const std::shared_ptr<const DiscreteFunctionVariant>& lambda,
                         const std::shared_ptr<const DiscreteFunctionVariant>& mu,
                         const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
                         const std::shared_ptr<const DiscreteFunctionVariant>& p_inf) const
  {
    std::shared_ptr sub_m = getCommonMesh({rho, u, E, CG});
    std::shared_ptr<const DiscreteFunctionVariant> sub_rho = rho;
    std::shared_ptr<const DiscreteFunctionVariant> sub_u   = u;
    std::shared_ptr<const DiscreteFunctionVariant> sub_E   = E;
    std::shared_ptr<const DiscreteFunctionVariant> sub_CG  = CG;

    std::shared_ptr<const DiscreteFunctionVariant> sub_lambda = lambda;
    std::shared_ptr<const DiscreteFunctionVariant> sub_mu     = mu;
    std::shared_ptr<const DiscreteFunctionVariant> sub_gamma  = gamma;
    std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf  = p_inf;

    double sum_dt = 0;
    while(sum_dt < dt){
      const auto& [sigma, p, aL, aT] = apply_state_law<MeshType>(law,sub_rho,sub_u,sub_E,sub_CG,sub_lambda,sub_mu,sub_gamma,sub_p_inf);

      const DiscreteScalarFunction& aL_d = aL->template get<DiscreteScalarFunction>();
      const DiscreteScalarFunction& aT_d = aT->template get<DiscreteScalarFunction>();
      const std::shared_ptr<const DiscreteFunctionVariant>& c =
        std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);

      double sub_dt = CFL * hyperelastic_dt(c);
      if(sum_dt + sub_dt > dt){
	      sub_dt = dt - sum_dt;
      }

      auto [sub_ur2, sub_Fjr2]       = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, p, bc_descriptor_list, CR_ur, CR_Fjr, map2);
      auto [sub_ur2_o1, sub_Fjr2_o1] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, bc_descriptor_list, CR_ur, CR_Fjr, map2);

      std::tie(sub_m, sub_rho, sub_u, sub_E, sub_CG) = apply_fluxes(sub_dt, sub_rho, sub_u, sub_E, sub_CG, sub_ur2, sub_ur2_o1, sub_Fjr2, sub_Fjr2_o1);

      sub_lambda = shallowCopy(sub_m,sub_lambda);
      sub_mu     = shallowCopy(sub_m,sub_mu);
      sub_gamma  = shallowCopy(sub_m,sub_gamma);
      sub_p_inf  = shallowCopy(sub_m,sub_p_inf);

      sum_dt += sub_dt;
    }

    return {sub_m, sub_rho, sub_u, sub_E, sub_CG};
  }

  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>>
  compute_sub_iterations(const SolverType& solver_type,
                         const size_t& law,
                         const double& dt,
                         const size_t& q,
                         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,
                         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                         NodeValue<Rd> CR_ur,
                         NodeValuePerCell<Rd> CR_Fjr,
                         std::vector<NodeId> map2,
                         const std::shared_ptr<const DiscreteFunctionVariant>& lambda,
                         const std::shared_ptr<const DiscreteFunctionVariant>& mu,
                         const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
                         const std::shared_ptr<const DiscreteFunctionVariant>& p_inf) const
  {
    std::shared_ptr sub_m = getCommonMesh({rho, u, E, CG});
    std::shared_ptr<const DiscreteFunctionVariant> sub_rho = rho;
    std::shared_ptr<const DiscreteFunctionVariant> sub_u   = u;
    std::shared_ptr<const DiscreteFunctionVariant> sub_E   = E;
    std::shared_ptr<const DiscreteFunctionVariant> sub_CG  = CG;

    std::shared_ptr<const DiscreteFunctionVariant> sub_lambda = lambda;
    std::shared_ptr<const DiscreteFunctionVariant> sub_mu     = mu;
    std::shared_ptr<const DiscreteFunctionVariant> sub_gamma  = gamma;
    std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf  = p_inf;

    double sub_dt = dt / q;
    double sum_dt = 0.;
    for(size_t i = 0; i < q; ++i){

      if( i == q - 1){
        sub_dt = dt - sum_dt;
      }
      
      const auto& [sigma, p, aL, aT] = apply_state_law<MeshType>(law,sub_rho,sub_u,sub_E,sub_CG,sub_lambda,sub_mu,sub_gamma,sub_p_inf);

      const DiscreteScalarFunction& aL_d = aL->template get<DiscreteScalarFunction>();
      const DiscreteScalarFunction& aT_d = aT->template get<DiscreteScalarFunction>();
      const std::shared_ptr<const DiscreteFunctionVariant>& c =
        std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);

      auto [sub_ur2, sub_Fjr2]       = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, p, bc_descriptor_list, CR_ur, CR_Fjr, map2);
      auto [sub_ur2_o1, sub_Fjr2_o1] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, bc_descriptor_list, CR_ur, CR_Fjr, map2);

      std::tie(sub_m, sub_rho, sub_u, sub_E, sub_CG) = apply_fluxes(sub_dt, sub_rho, sub_u, sub_E, sub_CG, sub_ur2, sub_ur2_o1, sub_Fjr2, sub_Fjr2_o1);

      sub_lambda = shallowCopy(sub_m,sub_lambda);
      sub_mu     = shallowCopy(sub_m,sub_mu);
      sub_gamma  = shallowCopy(sub_m,sub_gamma);
      sub_p_inf  = shallowCopy(sub_m,sub_p_inf);

      sum_dt += sub_dt;

    }

    return {sub_m, sub_rho, sub_u, sub_E, sub_CG};
  }

  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>>
  compute_sub_iterations_midpoint_method(const SolverType& solver_type,
                                         const size_t& law,
                                         const double& CFL,
                                         const double& dt,
                                         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,
                                         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                                         NodeValue<Rd> CR_ur,
                                         NodeValuePerCell<Rd> CR_Fjr,
                                         std::vector<NodeId> map2,
                                         const std::shared_ptr<const DiscreteFunctionVariant>& lambda,
                                         const std::shared_ptr<const DiscreteFunctionVariant>& mu,
                                         const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
                                         const std::shared_ptr<const DiscreteFunctionVariant>& p_inf) const
  {
    std::shared_ptr sub_m = getCommonMesh({rho, u});
    std::shared_ptr<const DiscreteFunctionVariant> sub_rho = rho;
    std::shared_ptr<const DiscreteFunctionVariant> sub_u   = u;
    std::shared_ptr<const DiscreteFunctionVariant> sub_E   = E;
    std::shared_ptr<const DiscreteFunctionVariant> sub_CG  = CG;

    double sum_dt = 0;
    while(sum_dt < dt){
      const std::shared_ptr<const DiscreteFunctionVariant> sub_lambda = shallowCopy(sub_m,lambda); 
      const std::shared_ptr<const DiscreteFunctionVariant> sub_mu     = shallowCopy(sub_m,mu); 
      const std::shared_ptr<const DiscreteFunctionVariant> sub_gamma  = shallowCopy(sub_m,gamma);
      const std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf  = shallowCopy(sub_m,p_inf);       
      const auto& [sigma, p, aL, aT] = apply_state_law<MeshType>(law,sub_rho,sub_u,sub_E,sub_CG,sub_lambda,sub_mu,sub_gamma,sub_p_inf);

      const DiscreteScalarFunction& aL_d = aL->template get<DiscreteScalarFunction>();
      const DiscreteScalarFunction& aT_d = aT->template get<DiscreteScalarFunction>();
      const std::shared_ptr<const DiscreteFunctionVariant>& c =
        std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);

      double sub_dt = CFL * hyperelastic_dt(c);
      if(sum_dt + sub_dt > dt){
	      sub_dt = dt - sum_dt;
      }

      auto [sub_ur, sub_Fjr]       = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, p, bc_descriptor_list, CR_ur, CR_Fjr, map2);
      auto [sub_ur_o1, sub_Fjr_o1] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, bc_descriptor_list, CR_ur, CR_Fjr, map2);

      auto [sub_m_app, sub_rho_app, sub_u_app, sub_E_app, sub_CG_app] = apply_fluxes(sub_dt/2., sub_rho, sub_u, sub_E, sub_CG, sub_ur, sub_ur_o1, sub_Fjr, sub_Fjr_o1);

      const std::shared_ptr<const DiscreteFunctionVariant> sub_lambda_app = shallowCopy(sub_m_app,lambda); 
      const std::shared_ptr<const DiscreteFunctionVariant> sub_mu_app     = shallowCopy(sub_m_app,mu); 
      const std::shared_ptr<const DiscreteFunctionVariant> sub_gamma_app  = shallowCopy(sub_m_app,gamma);
      const std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf_app  = shallowCopy(sub_m_app,p_inf);
      const auto& [sigma_app, p_app, aL_app, aT_app] = apply_state_law<MeshType>(law,sub_rho_app,sub_u_app,sub_E_app,sub_CG_app,sub_lambda_app,sub_mu_app,sub_gamma_app,sub_p_inf_app);

      auto [sub_ur_app, sub_Fjr_app]       = compute_fluxes(solver_type, sub_rho_app, aL_app, aT_app, sub_u_app, sigma_app, p_app, bc_descriptor_list, CR_ur, CR_Fjr, map2);
      auto [sub_ur_o1_app, sub_Fjr_o1_app] = compute_fluxes(solver_type, sub_rho_app, aL_app, aT_app, sub_u_app, sigma_app, bc_descriptor_list, CR_ur, CR_Fjr, map2);

      std::tie(sub_m, sub_rho, sub_u, sub_E, sub_CG) = order2_apply_fluxes(sub_dt, sub_rho, sub_u, sub_E, sub_CG, sub_rho_app, sub_CG_app, sub_ur_app, sub_ur_o1_app, sub_Fjr_app, sub_Fjr_o1_app); 

      sum_dt += sub_dt;
    }

    return {sub_m, sub_rho, sub_u, sub_E, sub_CG};
  }

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>>
compute_sub_iterations_midpoint_method(const SolverType& solver_type,
                              const size_t& law,
                              const double& dt,
                              const size_t q,
                              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,
                              const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                              NodeValue<Rd> CR_ur,
                              NodeValuePerCell<Rd> CR_Fjr,
                              std::vector<NodeId> map2,
                              const std::shared_ptr<const DiscreteFunctionVariant>& lambda,
                              const std::shared_ptr<const DiscreteFunctionVariant>& mu,
                              const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
                              const std::shared_ptr<const DiscreteFunctionVariant>& p_inf) const
{
  std::shared_ptr sub_m = getCommonMesh({rho, u});
  std::shared_ptr<const DiscreteFunctionVariant> sub_rho = rho;
  std::shared_ptr<const DiscreteFunctionVariant> sub_u   = u;
  std::shared_ptr<const DiscreteFunctionVariant> sub_E   = E;
  std::shared_ptr<const DiscreteFunctionVariant> sub_CG  = CG;

  double sub_dt = dt / q;
  double sum_dt = 0;
  for(size_t i = 0; i < q; ++i){

    if( i == q - 1){
      sub_dt = dt - sum_dt;
    }

    const std::shared_ptr<const DiscreteFunctionVariant> sub_lambda = shallowCopy(sub_m,lambda); 
    const std::shared_ptr<const DiscreteFunctionVariant> sub_mu     = shallowCopy(sub_m,mu); 
    const std::shared_ptr<const DiscreteFunctionVariant> sub_gamma  = shallowCopy(sub_m,gamma);
    const std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf  = shallowCopy(sub_m,p_inf);       
    const auto& [sigma, p, aL, aT] = apply_state_law<MeshType>(law,sub_rho,sub_u,sub_E,sub_CG,sub_lambda,sub_mu,sub_gamma,sub_p_inf);

    const DiscreteScalarFunction& aL_d = aL->template get<DiscreteScalarFunction>();
    const DiscreteScalarFunction& aT_d = aT->template get<DiscreteScalarFunction>();
    const std::shared_ptr<const DiscreteFunctionVariant>& c = std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);

    auto [sub_ur, sub_Fjr]       = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, p, bc_descriptor_list, CR_ur, CR_Fjr, map2);
    auto [sub_ur_o1, sub_Fjr_o1] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, bc_descriptor_list, CR_ur, CR_Fjr, map2);

    auto [sub_m_app, sub_rho_app, sub_u_app, sub_E_app, sub_CG_app] = apply_fluxes(sub_dt/2., sub_rho, sub_u, sub_E, sub_CG, sub_ur, sub_ur_o1, sub_Fjr, sub_Fjr_o1);

    const std::shared_ptr<const DiscreteFunctionVariant> sub_lambda_app = shallowCopy(sub_m_app,lambda); 
    const std::shared_ptr<const DiscreteFunctionVariant> sub_mu_app     = shallowCopy(sub_m_app,mu); 
    const std::shared_ptr<const DiscreteFunctionVariant> sub_gamma_app  = shallowCopy(sub_m_app,gamma);
    const std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf_app  = shallowCopy(sub_m_app,p_inf);
    const auto& [sigma_app, p_app, aL_app, aT_app] = apply_state_law<MeshType>(law,sub_rho_app,sub_u_app,sub_E_app,sub_CG_app,sub_lambda_app,sub_mu_app,sub_gamma_app,sub_p_inf_app);

    auto [sub_ur_app, sub_Fjr_app]       = compute_fluxes(solver_type, sub_rho_app, aL_app, aT_app, sub_u_app, sigma_app, p_app, bc_descriptor_list, CR_ur, CR_Fjr, map2);
    auto [sub_ur_o1_app, sub_Fjr_o1_app] = compute_fluxes(solver_type, sub_rho_app, aL_app, aT_app, sub_u_app, sigma_app, bc_descriptor_list, CR_ur, CR_Fjr, map2);

    std::tie(sub_m, sub_rho, sub_u, sub_E, sub_CG) = order2_apply_fluxes(sub_dt, sub_rho, sub_u, sub_E, sub_CG, sub_rho_app, sub_CG_app, sub_ur_app, sub_ur_o1_app, sub_Fjr_app, sub_Fjr_o1_app); 

    sum_dt += sub_dt;
  }

  return {sub_m, sub_rho, sub_u, sub_E, sub_CG};
}

  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>,
             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>>
  apply(const SolverType& solver_type,
        const size_t& law,
        const double& dt1,
        const size_t& q,
        const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
        const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
        const std::shared_ptr<const DiscreteFunctionVariant>& u1,
        const std::shared_ptr<const DiscreteFunctionVariant>& u2,
        const std::shared_ptr<const DiscreteFunctionVariant>& E1,
        const std::shared_ptr<const DiscreteFunctionVariant>& E2,
        const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
        const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
        const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
        const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
        const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
        const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
        const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
        const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
        const std::shared_ptr<const DiscreteFunctionVariant>& p1,
        const std::shared_ptr<const DiscreteFunctionVariant>& p2,
        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
        const std::shared_ptr<const DiscreteFunctionVariant>& lambda1,
        const std::shared_ptr<const DiscreteFunctionVariant>& mu1,
        const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
        const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
        const std::shared_ptr<const DiscreteFunctionVariant>& gamma1,
        const std::shared_ptr<const DiscreteFunctionVariant>& p_inf1,
        const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
        const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2) const
  {
    std::shared_ptr m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
    std::shared_ptr m1     = getCommonMesh({rho1, aL1, aT1, u1, sigma1});

    auto [map1, map2] = _computeMapping(m1,m2,bc_descriptor_list1,bc_descriptor_list2);

    auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr]                   = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,p1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,p2,bc_descriptor_list2,map1,map2);
    auto [ur1_o1, Fjr1_o1, ur2_o1, Fjr2_o1, CR_ur_o1, CR_Fjr_o1] = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,bc_descriptor_list2,map1,map2);
    
    auto [m1_app, rho1_app, u1_app, E1_app, CG1_app]   = apply_fluxes(dt1/2., rho1, u1, E1, CG1, ur1, ur1_o1, Fjr1, Fjr1_o1);

    const std::shared_ptr<const DiscreteFunctionVariant>& lambda1_app = shallowCopy(m1_app,lambda1);
    const std::shared_ptr<const DiscreteFunctionVariant>& mu1_app = shallowCopy(m1_app,mu1);
    const std::shared_ptr<const DiscreteFunctionVariant>& gamma1_app = shallowCopy(m1_app,gamma1);
    const std::shared_ptr<const DiscreteFunctionVariant>& p_inf1_app = shallowCopy(m1_app,p_inf1);
    const auto& [sigma1_app, p1_app, aL1_app, aT1_app] = apply_state_law<MeshType>(law, rho1_app, u1_app, E1_app, CG1_app, lambda1_app, mu1_app, gamma1_app, p_inf1_app);

    auto [m2_app, rho2_app, u2_app, E2_app, CG2_app]   = compute_sub_iterations(solver_type,law,dt1/2.,q,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur,CR_Fjr,map2,lambda2,mu2,gamma2,p_inf2);

    const std::shared_ptr<const DiscreteFunctionVariant>& lambda2_app = shallowCopy(m2_app,lambda2);
    const std::shared_ptr<const DiscreteFunctionVariant>& mu2_app = shallowCopy(m2_app,mu2);
    const std::shared_ptr<const DiscreteFunctionVariant>& gamma2_app = shallowCopy(m2_app,gamma2);
    const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2_app = shallowCopy(m2_app,p_inf2);
    const auto& [sigma2_app, p2_app, aL2_app, aT2_app] = apply_state_law<MeshType>(law, rho2_app, u2_app, E2_app, CG2_app, lambda2_app, mu2_app, gamma2_app, p_inf2_app);

    auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CR_ur_app, CR_Fjr_app]                   = compute_fluxes(solver_type,rho1_app,aL1_app,aT1_app,u1_app,sigma1_app,p1_app,bc_descriptor_list1,rho2_app,aL2_app,aT2_app,u2_app,sigma2_app,p2_app,bc_descriptor_list2,map1,map2);
    auto [ur1_o1_app, Fjr1_o1_app, ur2_o1_app, Fjr2_o1_app, CR_ur_o1_app, CR_Fjr_o1_app] = compute_fluxes(solver_type,rho1_app,aL1_app,aT1_app,u1_app,sigma1_app,bc_descriptor_list1,rho2_app,aL2_app,aT2_app,u2_app,sigma2_app,bc_descriptor_list2,map1,map2);

    const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = order2_apply_fluxes(dt1, rho1, u1, E1, CG1, rho1_app, CG1_app, ur1_app, ur1_o1_app, Fjr1_app, Fjr1_o1_app);
    auto [new_m2, new_rho2, new_u2, new_E2, new_CG2]       = compute_sub_iterations_midpoint_method(solver_type,law,dt1,q,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur_app,CR_Fjr_app,map2,lambda2,mu2,gamma2,p_inf2);

    std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) = mesh_correction<MeshType>(new_m1, new_m2, new_rho2, new_u2, new_E2, new_CG2, map1, map2);

    return {new_m1, new_rho1, new_u1, new_E1, new_CG1, new_m2, new_rho2, new_u2, new_E2, new_CG2};
  }

  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>,
             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>>
  apply(const SolverType& solver_type,
        const size_t& law, 
        const double& dt1,
        const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
        const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
        const std::shared_ptr<const DiscreteFunctionVariant>& u1,
        const std::shared_ptr<const DiscreteFunctionVariant>& u2,
        const std::shared_ptr<const DiscreteFunctionVariant>& E1,
        const std::shared_ptr<const DiscreteFunctionVariant>& E2,
        const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
        const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
        const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
        const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
        const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
        const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
        const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
        const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
        const std::shared_ptr<const DiscreteFunctionVariant>& p1,
        const std::shared_ptr<const DiscreteFunctionVariant>& p2,
        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
        const std::shared_ptr<const DiscreteFunctionVariant>& lambda1,
        const std::shared_ptr<const DiscreteFunctionVariant>& mu1,
        const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
        const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
        const std::shared_ptr<const DiscreteFunctionVariant>& gamma1,
        const std::shared_ptr<const DiscreteFunctionVariant>& p_inf1,
        const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
        const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2) const
  {
    std::shared_ptr m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
    std::shared_ptr m1     = getCommonMesh({rho1, aL1, aT1, u1, sigma1});

    auto [map1, map2] = _computeMapping(m1,m2,bc_descriptor_list1,bc_descriptor_list2);

    auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr]                   = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,p1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,p2,bc_descriptor_list2,map1,map2);
    auto [ur1_o1, Fjr1_o1, ur2_o1, Fjr2_o1, CR_ur_o1, CR_Fjr_o1] = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,bc_descriptor_list2,map1,map2);
    
    auto [m1_app, rho1_app, u1_app, E1_app, CG1_app]   = apply_fluxes(dt1/2., rho1, u1, E1, CG1, ur1, ur1_o1, Fjr1, Fjr1_o1);

    const std::shared_ptr<const DiscreteFunctionVariant>& lambda1_app = shallowCopy(m1_app,lambda1);
    const std::shared_ptr<const DiscreteFunctionVariant>& mu1_app = shallowCopy(m1_app,mu1);
    const std::shared_ptr<const DiscreteFunctionVariant>& gamma1_app = shallowCopy(m1_app,gamma1);
    const std::shared_ptr<const DiscreteFunctionVariant>& p_inf1_app = shallowCopy(m1_app,p_inf1);
    const auto& [sigma1_app, p1_app, aL1_app, aT1_app] = apply_state_law<MeshType>(law, rho1_app, u1_app, E1_app, CG1_app, lambda1_app, mu1_app, gamma1_app, p_inf1_app);

    auto [m2_app, rho2_app, u2_app, E2_app, CG2_app]   = compute_sub_iterations(solver_type,law,0.4,dt1/2.,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur,CR_Fjr,map2,lambda2,mu2,gamma2,p_inf2);

    const std::shared_ptr<const DiscreteFunctionVariant>& lambda2_app = shallowCopy(m2_app,lambda2);
    const std::shared_ptr<const DiscreteFunctionVariant>& mu2_app = shallowCopy(m2_app,mu2);
    const std::shared_ptr<const DiscreteFunctionVariant>& gamma2_app = shallowCopy(m2_app,gamma2);
    const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2_app = shallowCopy(m2_app,p_inf2);
    const auto& [sigma2_app, p2_app, aL2_app, aT2_app] = apply_state_law<MeshType>(law, rho2_app, u2_app, E2_app, CG2_app, lambda2_app, mu2_app, gamma2_app, p_inf2_app);

    auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CR_ur_app, CR_Fjr_app]                   = compute_fluxes(solver_type,rho1_app,aL1_app,aT1_app,u1_app,sigma1_app,p1_app,bc_descriptor_list1,rho2_app,aL2_app,aT2_app,u2_app,sigma2_app,p2_app,bc_descriptor_list2,map1,map2);
    auto [ur1_o1_app, Fjr1_o1_app, ur2_o1_app, Fjr2_o1_app, CR_ur_o1_app, CR_Fjr_o1_app] = compute_fluxes(solver_type,rho1_app,aL1_app,aT1_app,u1_app,sigma1_app,bc_descriptor_list1,rho2_app,aL2_app,aT2_app,u2_app,sigma2_app,bc_descriptor_list2,map1,map2);

    const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = order2_apply_fluxes(dt1, rho1, u1, E1, CG1, rho1_app, CG1_app, ur1_app, ur1_o1_app, Fjr1_app, Fjr1_o1_app);
    auto [new_m2, new_rho2, new_u2, new_E2, new_CG2]       = compute_sub_iterations_midpoint_method(solver_type,law,0.4,dt1,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur_app,CR_Fjr_app,map2,lambda2,mu2,gamma2,p_inf2);

    std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) = mesh_correction<MeshType>(new_m1, new_m2, new_rho2, new_u2, new_E2, new_CG2, map1, map2);

    return {new_m1, new_rho1, new_u1, new_E1, new_CG1, new_m2, new_rho2, new_u2, new_E2, new_CG2};
  }

  Order2LocalDtHyperelasticSolver()                            = default;
  Order2LocalDtHyperelasticSolver(Order2LocalDtHyperelasticSolver&&) = default;
  ~Order2LocalDtHyperelasticSolver()                           = default;
};

template <MeshConcept MeshType>
void
Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applyPressureBC(
  const BoundaryConditionList& bc_list,
  const MeshType& mesh,
  NodeValue<Rd>& br) const
{
  for (const auto& boundary_condition : bc_list) {
    std::visit(
      [&](auto&& bc) {
        using T = std::decay_t<decltype(bc)>;
        if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
          MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
          if constexpr (Dimension == 1) {
            const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();

            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

            const auto& node_list  = bc.nodeList();
            const auto& value_list = bc.valueList();
            parallel_for(
              node_list.size(), PUGS_LAMBDA(size_t i_node) {
                const NodeId node_id       = node_list[i_node];
                const auto& node_cell_list = node_to_cell_matrix[node_id];
                Assert(node_cell_list.size() == 1);

                CellId node_cell_id              = node_cell_list[0];
                size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);

                br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
              });
          } else {
            const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();

            const auto& face_to_cell_matrix               = mesh.connectivity().faceToCellMatrix();
            const auto& face_to_node_matrix               = mesh.connectivity().faceToNodeMatrix();
            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
            const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();

            const auto& face_list  = bc.faceList();
            const 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 auto& face_cell_list = face_to_cell_matrix[face_id];
              Assert(face_cell_list.size() == 1);

              CellId face_cell_id              = face_cell_list[0];
              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);

              const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;

              const auto& face_nodes = face_to_node_matrix[face_id];

              for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
                NodeId node_id = face_nodes[i_node];
                br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
              }
            }
          }
        }
      },
      boundary_condition);
  }
}

template <MeshConcept MeshType>
void
Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applyNormalStressBC(
  const BoundaryConditionList& bc_list,
  const MeshType& mesh,
  NodeValue<Rd>& br) const
{
  for (const auto& boundary_condition : bc_list) {
    std::visit(
      [&](auto&& bc) {
        using T = std::decay_t<decltype(bc)>;
        if constexpr (std::is_same_v<NormalStressBoundaryCondition, T>) {
          MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
          if constexpr (Dimension == 1) {
            const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();

            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

            const auto& node_list  = bc.nodeList();
            const auto& value_list = bc.valueList();
            parallel_for(
              node_list.size(), PUGS_LAMBDA(size_t i_node) {
                const NodeId node_id       = node_list[i_node];
                const auto& node_cell_list = node_to_cell_matrix[node_id];
                Assert(node_cell_list.size() == 1);

                CellId node_cell_id              = node_cell_list[0];
                size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);

                br[node_id] += value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
              });
          } else {
            const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();

            const auto& face_to_cell_matrix               = mesh.connectivity().faceToCellMatrix();
            const auto& face_to_node_matrix               = mesh.connectivity().faceToNodeMatrix();
            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
            const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();

            const auto& face_list  = bc.faceList();
            const 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 auto& face_cell_list = face_to_cell_matrix[face_id];
              Assert(face_cell_list.size() == 1);

              CellId face_cell_id              = face_cell_list[0];
              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);

              const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;

              const auto& face_nodes = face_to_node_matrix[face_id];

              for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
                NodeId node_id = face_nodes[i_node];
                br[node_id] += sign * value_list[i_face] * Nlr(face_id, i_node);
              }
            }
          }
        }
      },
      boundary_condition);
  }
}

template <MeshConcept MeshType>
void
Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applySymmetryBC(
  const BoundaryConditionList& bc_list,
  NodeValue<Rdxd>& Ar,
  NodeValue<Rd>& br) const
{
  for (const auto& boundary_condition : bc_list) {
    std::visit(
      [&](auto&& bc) {
        using T = std::decay_t<decltype(bc)>;
        if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
          const Rd& n = bc.outgoingNormal();

          const Rdxd I   = identity;
          const Rdxd nxn = tensorProduct(n, n);
          const Rdxd P   = I - nxn;

          const Array<const NodeId>& node_list = bc.nodeList();
          parallel_for(
            bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
              const NodeId r = node_list[r_number];

              Ar[r] = P * Ar[r] * P + nxn;
              br[r] = P * br[r];
            });
        }
      },
      boundary_condition);
  }
}

template <MeshConcept MeshType>
void
Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applyVelocityBC(
  const BoundaryConditionList& bc_list,
  NodeValue<Rdxd>& Ar,
  NodeValue<Rd>& br) const
{
  for (const auto& boundary_condition : bc_list) {
    std::visit(
      [&](auto&& bc) {
        using T = std::decay_t<decltype(bc)>;
        if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
          const auto& node_list  = bc.nodeList();
          const auto& value_list = bc.valueList();

          parallel_for(
            node_list.size(), PUGS_LAMBDA(size_t i_node) {
              NodeId node_id    = node_list[i_node];
              const auto& value = value_list[i_node];

              Ar[node_id] = identity;
              br[node_id] = value;
            });
        } else if constexpr (std::is_same_v<FixedBoundaryCondition, T>) {
          const auto& node_list = bc.nodeList();
          parallel_for(
            node_list.size(), PUGS_LAMBDA(size_t i_node) {
              NodeId node_id = node_list[i_node];

              Ar[node_id] = identity;
              br[node_id] = zero;
            });
        }
      },
      boundary_condition);
  }
}

template <MeshConcept MeshType>
void
Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC(
  NodeValue<Rdxd>& Ar1,
  NodeValue<Rdxd>& Ar2,
  NodeValue<Rd>& br1,
  NodeValue<Rd>& br2,
  const std::vector<NodeId>& map1,
  const std::vector<NodeId>& map2) const
{
  size_t n = map1.size();

  for (size_t i = 0; i < n; i++) {
    NodeId node_id1 = map1[i];
    NodeId node_id2 = map2[i];

    Ar1[node_id1] += Ar2[node_id2];
    Ar2[node_id2] = Ar1[node_id1];

    br1[node_id1] += br2[node_id2];
    br2[node_id2] = br1[node_id1];
  }
}

template <MeshConcept MeshType>
void 
Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC2(NodeValue<Rdxd>& Ar1,
										      NodeValue<Rdxd>& Ar2,
										      NodeValue<Rd>& ur1,
										      NodeValue<Rd>& ur2,
										      const std::vector<NodeId>& map1,
										      const std::vector<NodeId>& map2) const
{
  size_t n = map1.size();
  Rd lambda;

  for(size_t i = 0; i<n; i++){

    NodeId node_id1 = map1[i];
    NodeId node_id2 = map2[i];
    
    lambda = inverse(inverse(Ar2[node_id2]) + inverse(Ar1[node_id1]))*(ur2[node_id2]-ur1[node_id1]); 

    ur1[node_id1] = ur1[node_id1] + inverse(Ar1[node_id1])*lambda;
    ur2[node_id2] = ur2[node_id2] - inverse(Ar2[node_id2])*lambda; 
  }
}

template <MeshConcept MeshType>
void
Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC(
  const MeshType& mesh,
  NodeValue<Rd>& ur,
  NodeValue<Rd>& CR_ur,
  NodeValuePerCell<Rd>& Fjr,
  NodeValuePerCell<Rd>& CR_Fjr,
  const std::vector<NodeId>& map2) const
{
  const size_t& n = map2.size();

  const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
  const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

  for (size_t i = 0; i < n; i++) {
    const NodeId node_id                       = map2[i];
    const auto& node_to_cell                   = node_to_cell_matrix[node_id];
    const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);

    for (size_t j = 0; j < node_to_cell.size(); j++) {
      const CellId J       = node_to_cell[j];
      const unsigned int R = node_local_number_in_its_cells[j];
      Fjr(J, R)            = CR_Fjr(J, R);
    }
    ur[node_id] = CR_ur[node_id];
  }
}

template <MeshConcept MeshType>
class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::FixedBoundaryCondition
{
 private:
  const MeshNodeBoundary m_mesh_node_boundary;

 public:
  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  FixedBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}

  ~FixedBoundaryCondition() = default;
};

template <MeshConcept MeshType>
class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::PressureBoundaryCondition
{
 private:
  const MeshFaceBoundary m_mesh_face_boundary;
  const Array<const double> m_value_list;

 public:
  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_face_boundary.faceList();
  }

  const Array<const double>&
  valueList() const
  {
    return m_value_list;
  }

  PressureBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const double>& value_list)
    : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
  {}

  ~PressureBoundaryCondition() = default;
};

template <>
class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<Mesh<1>>::PressureBoundaryCondition
{
 private:
  const MeshNodeBoundary m_mesh_node_boundary;
  const Array<const double> m_value_list;

 public:
  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const double>&
  valueList() const
  {
    return m_value_list;
  }

  PressureBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const double>& value_list)
    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
  {}

  ~PressureBoundaryCondition() = default;
};

template <MeshConcept MeshType>
class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::NormalStressBoundaryCondition
{
 private:
  const MeshFaceBoundary m_mesh_face_boundary;
  const Array<const Rdxd> m_value_list;

 public:
  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_face_boundary.faceList();
  }

  const Array<const Rdxd>&
  valueList() const
  {
    return m_value_list;
  }

  NormalStressBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const Rdxd>& value_list)
    : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
  {}

  ~NormalStressBoundaryCondition() = default;
};

template <>
class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<Mesh<1>>::NormalStressBoundaryCondition
{
 private:
  const MeshNodeBoundary m_mesh_node_boundary;
  const Array<const Rdxd> m_value_list;

 public:
  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const Rdxd>&
  valueList() const
  {
    return m_value_list;
  }

  NormalStressBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const Rdxd>& value_list)
    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
  {}

  ~NormalStressBoundaryCondition() = default;
};

template <MeshConcept MeshType>
class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::VelocityBoundaryCondition
{
 private:
  const MeshNodeBoundary m_mesh_node_boundary;

  const Array<const TinyVector<Dimension>> m_value_list;

 public:
  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const TinyVector<Dimension>>&
  valueList() const
  {
    return m_value_list;
  }

  VelocityBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
                            const Array<const TinyVector<Dimension>>& value_list)
    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
  {}

  ~VelocityBoundaryCondition() = default;
};

template <MeshConcept MeshType>
class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::SymmetryBoundaryCondition
{
 public:
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;

 public:
  const Rd&
  outgoingNormal() const
  {
    return m_mesh_flat_node_boundary.outgoingNormal();
  }

  size_t
  numberOfNodes() const
  {
    return m_mesh_flat_node_boundary.nodeList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_flat_node_boundary.nodeList();
  }

  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary)
    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
  {
    ;
  }

  ~SymmetryBoundaryCondition() = default;
};

template <MeshConcept MeshType>
class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::CouplingBoundaryCondition
{
 private:
  const MeshNodeBoundary m_mesh_node_boundary;
  const size_t m_label;

 public:
  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  size_t
  label() const
  {
    return m_label;
  }

  CouplingBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const size_t label)
    : m_mesh_node_boundary{mesh_node_boundary}, m_label{label}
  {
    ;
  }

  ~CouplingBoundaryCondition() = default;
};

Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolverHandler(const std::shared_ptr<const MeshVariant>& i_mesh1,
                                                                   const std::shared_ptr<const MeshVariant>& i_mesh2)
{
  if (not i_mesh1) {
    throw NormalError("mesh1 not defined");
  }

  if (not i_mesh2) {
    throw NormalError("mesh2 not defined");
  }
  std::visit(
    [&](auto&& first_mesh, auto&& second_mesh) {
      using FirstMeshType  = mesh_type_t<decltype(first_mesh)>;
      using SecondMeshType = mesh_type_t<decltype(second_mesh)>;

      if constexpr (!std::is_same_v<FirstMeshType,
                                    SecondMeshType>) {   // <- ICI POUR LE TEST SUR LES TYPES DE MAILLAGES
        throw NormalError("incompatible mesh types");
      }
    },
    i_mesh1->variant(), i_mesh2->variant());

  std::visit(
    [&](auto&& mesh) {
      using MeshType = mesh_type_t<decltype(mesh)>;
      if constexpr (is_polygonal_mesh_v<MeshType>) {
        m_hyperelastic_solver = std::make_unique<Order2LocalDtHyperelasticSolver<MeshType>>();
      } else {
        throw NormalError("unexpected mesh type");
      }
    },
    i_mesh1->variant());
}