#include <language/modules/SchemeModule.hpp>

#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <language/modules/BinaryOperatorRegisterForVh.hpp>
#include <language/modules/MathFunctionRegisterForVh.hpp>
#include <language/modules/UnaryOperatorRegisterForVh.hpp>
#include <language/utils/BinaryOperatorProcessorBuilder.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
#include <mesh/IZoneDescriptor.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshRandomizer.hpp>
#include <mesh/MeshSmoother.hpp>
#include <scheme/AcousticSolver.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionDescriptorP0.hpp>
#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
#include <scheme/DiscreteFunctionIntegrator.hpp>
#include <scheme/DiscreteFunctionInterpoler.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/DiscreteFunctionVectorIntegrator.hpp>
#include <scheme/DiscreteFunctionVectorInterpoler.hpp>
#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
#include <scheme/FluxingAdvectionSolver.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
#include <scheme/FreeBoundaryConditionDescriptor.hpp>
#include <scheme/HyperelasticSolver.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <utils/Socket.hpp>

#include <memory>

SchemeModule::SchemeModule()
{
  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const DiscreteFunctionVariant>>);
  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IDiscreteFunctionDescriptor>>);
  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IQuadratureDescriptor>>);

  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>>);

  this->_addBuiltinFunction("P0", std::function(

                                    []() -> std::shared_ptr<const IDiscreteFunctionDescriptor> {
                                      return std::make_shared<DiscreteFunctionDescriptorP0>();
                                    }

                                    ));

  this->_addBuiltinFunction("P0Vector", std::function(

                                          []() -> std::shared_ptr<const IDiscreteFunctionDescriptor> {
                                            return std::make_shared<DiscreteFunctionDescriptorP0Vector>();
                                          }

                                          ));

  this->_addBuiltinFunction("Gauss", std::function(

                                       [](uint64_t degree) -> std::shared_ptr<const IQuadratureDescriptor> {
                                         return std::make_shared<GaussQuadratureDescriptor>(degree);
                                       }

                                       ));

  this->_addBuiltinFunction("GaussLobatto", std::function(

                                              [](uint64_t degree) -> std::shared_ptr<const IQuadratureDescriptor> {
                                                return std::make_shared<GaussLobattoQuadratureDescriptor>(degree);
                                              }

                                              ));

  this->_addBuiltinFunction("GaussLegendre", std::function(

                                               [](uint64_t degree) -> std::shared_ptr<const IQuadratureDescriptor> {
                                                 return std::make_shared<GaussLegendreQuadratureDescriptor>(degree);
                                               }

                                               ));

  this->_addBuiltinFunction("integrate",
                            std::function(

                              [](std::shared_ptr<const IMesh> mesh,
                                 const std::vector<std::shared_ptr<const IZoneDescriptor>>& integration_zone_list,
                                 std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor,
                                 std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
                                 const std::vector<FunctionSymbolId>& function_id_list)
                                -> std::shared_ptr<const DiscreteFunctionVariant> {
                                return std::make_shared<DiscreteFunctionVariant>(
                                  DiscreteFunctionVectorIntegrator{mesh, integration_zone_list, quadrature_descriptor,
                                                                   discrete_function_descriptor, function_id_list}
                                    .integrate());
                              }

                              ));

  this->_addBuiltinFunction("integrate",
                            std::function(

                              [](std::shared_ptr<const IMesh> mesh,
                                 std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor,
                                 std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
                                 const std::vector<FunctionSymbolId>& function_id_list)
                                -> std::shared_ptr<const DiscreteFunctionVariant> {
                                return std::make_shared<DiscreteFunctionVariant>(
                                  DiscreteFunctionVectorIntegrator{mesh, quadrature_descriptor,
                                                                   discrete_function_descriptor, function_id_list}
                                    .integrate());
                              }

                              ));

  this->_addBuiltinFunction("integrate",
                            std::function(

                              [](std::shared_ptr<const IMesh> mesh,
                                 const std::vector<std::shared_ptr<const IZoneDescriptor>>& integration_zone_list,
                                 std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor,
                                 const FunctionSymbolId& function_id)
                                -> std::shared_ptr<const DiscreteFunctionVariant> {
                                return std::make_shared<DiscreteFunctionVariant>(
                                  DiscreteFunctionIntegrator{mesh, integration_zone_list, quadrature_descriptor,
                                                             function_id}
                                    .integrate());
                              }

                              ));

  this->_addBuiltinFunction("integrate",
                            std::function(

                              [](std::shared_ptr<const IMesh> mesh,
                                 std::shared_ptr<const IQuadratureDescriptor> quadrature_descriptor,
                                 const FunctionSymbolId& function_id)
                                -> std::shared_ptr<const DiscreteFunctionVariant> {
                                return std::make_shared<DiscreteFunctionVariant>(
                                  DiscreteFunctionIntegrator{mesh, quadrature_descriptor, function_id}.integrate());
                              }

                              ));

  this->_addBuiltinFunction("interpolate",
                            std::function(

                              [](std::shared_ptr<const IMesh> mesh,
                                 const std::vector<std::shared_ptr<const IZoneDescriptor>>& interpolation_zone_list,
                                 std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
                                 const std::vector<FunctionSymbolId>& function_id_list)
                                -> std::shared_ptr<const DiscreteFunctionVariant> {
                                switch (discrete_function_descriptor->type()) {
                                case DiscreteFunctionType::P0: {
                                  if (function_id_list.size() != 1) {
                                    throw NormalError("invalid function descriptor type");
                                  }
                                  return std::make_shared<DiscreteFunctionVariant>(
                                    DiscreteFunctionInterpoler{mesh, interpolation_zone_list,
                                                               discrete_function_descriptor, function_id_list[0]}
                                      .interpolate());
                                }
                                case DiscreteFunctionType::P0Vector: {
                                  return std::make_shared<DiscreteFunctionVariant>(
                                    DiscreteFunctionVectorInterpoler{mesh, interpolation_zone_list,
                                                                     discrete_function_descriptor, function_id_list}
                                      .interpolate());
                                }
                                default: {
                                  throw NormalError("invalid function descriptor type");
                                }
                                }
                              }

                              ));

  this->_addBuiltinFunction("interpolate",
                            std::function(

                              [](std::shared_ptr<const IMesh> mesh,
                                 std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
                                 const std::vector<FunctionSymbolId>& function_id_list)
                                -> std::shared_ptr<const DiscreteFunctionVariant> {
                                switch (discrete_function_descriptor->type()) {
                                case DiscreteFunctionType::P0: {
                                  if (function_id_list.size() != 1) {
                                    throw NormalError("invalid function descriptor type");
                                  }
                                  return std::make_shared<DiscreteFunctionVariant>(
                                    DiscreteFunctionInterpoler{mesh, discrete_function_descriptor, function_id_list[0]}
                                      .interpolate());
                                }
                                case DiscreteFunctionType::P0Vector: {
                                  return std::make_shared<DiscreteFunctionVariant>(
                                    DiscreteFunctionVectorInterpoler{mesh, discrete_function_descriptor,
                                                                     function_id_list}
                                      .interpolate());
                                }
                                default: {
                                  throw NormalError("invalid function descriptor type");
                                }
                                }
                              }

                              ));

  this->_addBuiltinFunction("randomizeMesh",
                            std::function(

                              [](std::shared_ptr<const IMesh> p_mesh,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list) -> std::shared_ptr<const IMesh> {
                                MeshRandomizerHandler handler;
                                return handler.getRandomizedMesh(*p_mesh, bc_descriptor_list);
                              }

                              ));

  this->_addBuiltinFunction("randomizeMesh",
                            std::function(

                              [](std::shared_ptr<const IMesh> p_mesh,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const FunctionSymbolId& function_symbol_id) -> std::shared_ptr<const IMesh> {
                                MeshRandomizerHandler handler;
                                return handler.getRandomizedMesh(*p_mesh, bc_descriptor_list, function_symbol_id);
                              }

                              ));

  this->_addBuiltinFunction("smoothMesh", std::function(

                                            [](std::shared_ptr<const IMesh> p_mesh,
                                               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                                 bc_descriptor_list) -> std::shared_ptr<const IMesh> {
                                              MeshSmootherHandler handler;
                                              return handler.getSmoothedMesh(*p_mesh, bc_descriptor_list);
                                            }

                                            ));

  this->_addBuiltinFunction("smoothMesh",
                            std::function(

                              [](std::shared_ptr<const IMesh> p_mesh,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const FunctionSymbolId& function_symbol_id) -> std::shared_ptr<const IMesh> {
                                MeshSmootherHandler handler;
                                return handler.getSmoothedMesh(*p_mesh, bc_descriptor_list, function_symbol_id);
                              }

                              ));

  this->_addBuiltinFunction("fixed", std::function(

                                       [](std::shared_ptr<const IBoundaryDescriptor> boundary)
                                         -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                         return std::make_shared<FixedBoundaryConditionDescriptor>(boundary);
                                       }

                                       ));

  this->_addBuiltinFunction("axis", std::function(

                                      [](std::shared_ptr<const IBoundaryDescriptor> boundary)
                                        -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                        return std::make_shared<AxisBoundaryConditionDescriptor>(boundary);
                                      }

                                      ));

  this->_addBuiltinFunction("symmetry", std::function(

                                          [](std::shared_ptr<const IBoundaryDescriptor> boundary)
                                            -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                            return std::make_shared<SymmetryBoundaryConditionDescriptor>(boundary);
                                          }

                                          ));

  this->_addBuiltinFunction("pressure", std::function(

                                          [](std::shared_ptr<const IBoundaryDescriptor> boundary,
                                             const FunctionSymbolId& pressure_id)
                                            -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                            return std::make_shared<DirichletBoundaryConditionDescriptor>("pressure",
                                                                                                          boundary,
                                                                                                          pressure_id);
                                          }

                                          ));

  this->_addBuiltinFunction("normalstress",
                            std::function(

                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
                                 const FunctionSymbolId& normal_stress_id)
                                -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                return std::make_shared<DirichletBoundaryConditionDescriptor>("normal-stress", boundary,
                                                                                              normal_stress_id);
                              }

                              ));

  this->_addBuiltinFunction("velocity", std::function(

                                          [](std::shared_ptr<const IBoundaryDescriptor> boundary,
                                             const FunctionSymbolId& velocity_id)
                                            -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                            return std::make_shared<DirichletBoundaryConditionDescriptor>("velocity",
                                                                                                          boundary,
                                                                                                          velocity_id);
                                          }

                                          ));

  this->_addBuiltinFunction("external_fsi_velocity",
                            std::function(

                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
                                 const std::shared_ptr<const Socket>& socket)
                                -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                return std::make_shared<ExternalBoundaryConditionDescriptor>("external_fsi_velocity",
                                                                                             boundary, socket);
                              }

                              ));

  this->_addBuiltinFunction("glace_fluxes", std::function(

                                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
                                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
                                                 const std::shared_ptr<const DiscreteFunctionVariant>& c,
                                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                                   bc_descriptor_list)
                                                -> std::tuple<std::shared_ptr<const ItemValueVariant>,
                                                              std::shared_ptr<const SubItemValuePerItemVariant>> {
                                                return AcousticSolverHandler{getCommonMesh({rho, c, u, p})}
                                                  .solver()
                                                  .compute_fluxes(AcousticSolverHandler::SolverType::Glace, rho, c, u,
                                                                  p, bc_descriptor_list);
                                              }

                                              ));

  this->_addBuiltinFunction("glace_solver",
                            std::function(

                              [](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>& c,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
                                return AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
                                  .solver()
                                  .apply(AcousticSolverHandler::SolverType::Glace, dt, rho, u, E, c, p,
                                         bc_descriptor_list);
                              }

                              ));

  this->_addBuiltinFunction("eucclhyd_fluxes",
                            std::function(

                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& c,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list)
                                -> std::tuple<std::shared_ptr<const ItemValueVariant>,
                                              std::shared_ptr<const SubItemValuePerItemVariant>> {
                                return AcousticSolverHandler{getCommonMesh({rho, c, u, p})}
                                  .solver()
                                  .compute_fluxes(AcousticSolverHandler::SolverType::Eucclhyd, rho, c, u, p,
                                                  bc_descriptor_list);
                              }

                              ));

  this->_addBuiltinFunction("eucclhyd_solver",
                            std::function(

                              [](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>& c,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
                                return AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
                                  .solver()
                                  .apply(AcousticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, c, p,
                                         bc_descriptor_list);
                              }

                              ));

  this->_addBuiltinFunction("apply_acoustic_fluxes",
                            std::function(

                              [](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 ItemValueVariant>& ur,              //
                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,   //
                                 const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
                                return AcousticSolverHandler{getCommonMesh({rho, u, E})}   //
                                  .solver()
                                  .apply_fluxes(dt, rho, u, E, ur, Fjr);
                              }

                              ));

  this->_addBuiltinFunction("hyperelastic_eucclhyd_fluxes",
                            std::function(

                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list)
                                -> std::tuple<std::shared_ptr<const ItemValueVariant>,
                                              std::shared_ptr<const SubItemValuePerItemVariant>> {
                                return HyperelasticSolverHandler{getCommonMesh({rho, aL, aT, u, sigma})}
                                  .solver()
                                  .compute_fluxes(HyperelasticSolverHandler::SolverType::Eucclhyd, rho, aL, aT, u,
                                                  sigma, bc_descriptor_list);
                              }

                              ));

  this->_addBuiltinFunction("hyperelastic_eucclhyd_solver",
                            std::function(

                              [](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>& aL,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
                                return HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
                                  .solver()
                                  .apply(HyperelasticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, CG, aL, aT,
                                         sigma, bc_descriptor_list);
                              }

                              ));

  this->_addBuiltinFunction("hyperelastic_glace_fluxes",
                            std::function(

                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list)
                                -> std::tuple<std::shared_ptr<const ItemValueVariant>,
                                              std::shared_ptr<const SubItemValuePerItemVariant>> {
                                return HyperelasticSolverHandler{getCommonMesh({rho, aL, aT, u, sigma})}
                                  .solver()
                                  .compute_fluxes(HyperelasticSolverHandler::SolverType::Glace,   //
                                                  rho, aL, aT, u, sigma, bc_descriptor_list);
                              }

                              ));

  this->_addBuiltinFunction("hyperelastic_glace_solver",
                            std::function(

                              [](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>& aL,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
                                return HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
                                  .solver()
                                  .apply(HyperelasticSolverHandler::SolverType::Glace, dt, rho, u, E, CG, aL, aT, sigma,
                                         bc_descriptor_list);
                              }

                              ));

  this->_addBuiltinFunction("apply_hyperelastic_fluxes",
                            std::function(

                              [](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 SubItemValuePerItemVariant>& Fjr,   //
                                 const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
                                return HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG})}   //
                                  .solver()
                                  .apply_fluxes(dt, rho, u, E, CG, ur, Fjr);
                              }

                              ));

  this->_addBuiltinFunction("lagrangian",
                            std::function(

                              [](const std::shared_ptr<const IMesh>& mesh,
                                 const std::shared_ptr<const DiscreteFunctionVariant>& v)
                                -> std::shared_ptr<const DiscreteFunctionVariant> { return shallowCopy(mesh, v); }

                              ));

  this->_addBuiltinFunction("acoustic_dt", std::function(

                                             [](const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
                                               return acoustic_dt(c);
                                             }

                                             ));

  this->_addBuiltinFunction("cell_volume",
                            std::function(

                              [](const std::shared_ptr<const IMesh>& i_mesh)
                                -> std::shared_ptr<const DiscreteFunctionVariant> {
                                switch (i_mesh->dimension()) {
                                case 1: {
                                  constexpr size_t Dimension = 1;
                                  using MeshType             = Mesh<Connectivity<Dimension>>;
                                  std::shared_ptr<const MeshType> mesh =
                                    std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(i_mesh);

                                  return std::make_shared<DiscreteFunctionVariant>(
                                    DiscreteFunctionP0(mesh, MeshDataManager::instance().getMeshData(*mesh).Vj()));
                                }
                                case 2: {
                                  constexpr size_t Dimension = 2;
                                  using MeshType             = Mesh<Connectivity<Dimension>>;
                                  std::shared_ptr<const MeshType> mesh =
                                    std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(i_mesh);

                                  return std::make_shared<DiscreteFunctionVariant>(
                                    DiscreteFunctionP0(mesh, MeshDataManager::instance().getMeshData(*mesh).Vj()));
                                }
                                case 3: {
                                  constexpr size_t Dimension = 3;
                                  using MeshType             = Mesh<Connectivity<Dimension>>;
                                  std::shared_ptr<const MeshType> mesh =
                                    std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(i_mesh);

                                  return std::make_shared<DiscreteFunctionVariant>(
                                    DiscreteFunctionP0(mesh, MeshDataManager::instance().getMeshData(*mesh).Vj()));
                                }
                                default: {
                                  throw UnexpectedError("invalid mesh dimension");
                                }
                                }
                              }

                              ));

  this->_addBuiltinFunction("hyperelastic_dt", std::function(

                                                 [](const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
                                                   return hyperelastic_dt(c);
                                                 }

                                                 ));
  this->_addBuiltinFunction("fluxing_advection",
                            std::function(

                              [](const std::shared_ptr<const IMesh> new_mesh,
                                 const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& remapped_variables)
                                -> std::vector<std::shared_ptr<const DiscreteFunctionVariant>> {
                                return advectByFluxing(new_mesh, remapped_variables);
                              }

                              ));

  MathFunctionRegisterForVh{*this};
}

void
SchemeModule::registerOperators() const
{
  BinaryOperatorRegisterForVh{};
  UnaryOperatorRegisterForVh{};
}