#include <language/modules/SchemeModule.hpp>

#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussLobattoQuadratureDescriptor.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <language/algorithms/ElasticityDiamondAlgorithm.hpp>
#include <language/algorithms/Heat5PointsAlgorithm.hpp>
#include <language/algorithms/HeatDiamondAlgorithm.hpp>
#include <language/algorithms/HeatDiamondAlgorithm2.hpp>
#include <language/algorithms/ParabolicHeat.hpp>
#include <language/algorithms/UnsteadyElasticity.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 <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/DiscreteFunctionVectorIntegrator.hpp>
#include <scheme/DiscreteFunctionVectorInterpoler.hpp>
#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
#include <scheme/FreeBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunction.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/ScalarDiamondScheme.hpp>
#include <scheme/ScalarNodalScheme.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/VectorDiamondScheme.hpp>
#include <utils/Socket.hpp>

#include <memory>

SchemeModule::SchemeModule()
{
  this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IDiscreteFunction>>);
  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::make_shared<
                                    BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunctionDescriptor>()>>(
                                    []() -> std::shared_ptr<const IDiscreteFunctionDescriptor> {
                                      return std::make_shared<DiscreteFunctionDescriptorP0>();
                                    }

                                    ));

  this->_addBuiltinFunction("P0Vector",
                            std::make_shared<
                              BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunctionDescriptor>()>>(
                              []() -> std::shared_ptr<const IDiscreteFunctionDescriptor> {
                                return std::make_shared<DiscreteFunctionDescriptorP0Vector>();
                              }

                              ));

  this->_addBuiltinFunction("Gauss", std::make_shared<
                                       BuiltinFunctionEmbedder<std::shared_ptr<const IQuadratureDescriptor>(uint64_t)>>(
                                       [](uint64_t degree) -> std::shared_ptr<const IQuadratureDescriptor> {
                                         return std::make_shared<GaussQuadratureDescriptor>(degree);
                                       }

                                       ));

  this->_addBuiltinFunction("GaussLobatto",
                            std::make_shared<
                              BuiltinFunctionEmbedder<std::shared_ptr<const IQuadratureDescriptor>(uint64_t)>>(
                              [](uint64_t degree) -> std::shared_ptr<const IQuadratureDescriptor> {
                                return std::make_shared<GaussLobattoQuadratureDescriptor>(degree);
                              }

                              ));

  this->_addBuiltinFunction("GaussLegendre",
                            std::make_shared<
                              BuiltinFunctionEmbedder<std::shared_ptr<const IQuadratureDescriptor>(uint64_t)>>(
                              [](uint64_t degree) -> std::shared_ptr<const IQuadratureDescriptor> {
                                return std::make_shared<GaussLegendreQuadratureDescriptor>(degree);
                              }

                              ));

  this->_addBuiltinFunction("integrate",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IDiscreteFunction>(std::shared_ptr<const IMesh>,
                                                       const std::vector<std::shared_ptr<const IZoneDescriptor>>&,
                                                       std::shared_ptr<const IQuadratureDescriptor>,
                                                       std::shared_ptr<const IDiscreteFunctionDescriptor>,
                                                       const std::vector<FunctionSymbolId>&)>>(
                              [](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 IDiscreteFunction> {
                                return DiscreteFunctionVectorIntegrator{mesh, integration_zone_list,
                                                                        quadrature_descriptor,
                                                                        discrete_function_descriptor, function_id_list}
                                  .integrate();
                              }

                              ));

  this
    ->_addBuiltinFunction("integrate",
                          std::make_shared<BuiltinFunctionEmbedder<
                            std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IMesh>,
                                                                     std::shared_ptr<const IQuadratureDescriptor>,
                                                                     std::shared_ptr<const IDiscreteFunctionDescriptor>,
                                                                     const std::vector<FunctionSymbolId>&)>>(
                            [](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 IDiscreteFunction> {
                              return DiscreteFunctionVectorIntegrator{mesh, quadrature_descriptor,
                                                                      discrete_function_descriptor, function_id_list}
                                .integrate();
                            }

                            ));

  this->_addBuiltinFunction(
    "integrate",
    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
      const IDiscreteFunction>(std::shared_ptr<const IMesh>, const std::vector<std::shared_ptr<const IZoneDescriptor>>&,
                               std::shared_ptr<const IQuadratureDescriptor>, const FunctionSymbolId&)>>(
      [](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 IDiscreteFunction> {
        return DiscreteFunctionIntegrator{mesh, integration_zone_list, quadrature_descriptor, function_id}.integrate();
      }

      ));

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

                              ));

  this->_addBuiltinFunction(
    "interpolate",
    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
      const IDiscreteFunction>(std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunctionDescriptor>,
                               const std::vector<FunctionSymbolId>&)>>(
      [](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 IDiscreteFunction> {
        return DiscreteFunctionVectorInterpoler{mesh, discrete_function_descriptor, function_id_list}.interpolate();
      }

      ));

  this->_addBuiltinFunction("interpolate",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IDiscreteFunction>(std::shared_ptr<const IMesh>,
                                                       const std::vector<std::shared_ptr<const IZoneDescriptor>>&,
                                                       std::shared_ptr<const IDiscreteFunctionDescriptor>,
                                                       const std::vector<FunctionSymbolId>&)>>(
                              [](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 IDiscreteFunction> {
                                return DiscreteFunctionVectorInterpoler{mesh, interpolation_zone_list,
                                                                        discrete_function_descriptor, function_id_list}
                                  .interpolate();
                              }

                              ));

  this->_addBuiltinFunction(
    "interpolate",
    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
      const IDiscreteFunction>(std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunctionDescriptor>,
                               const FunctionSymbolId&)>>(
      [](std::shared_ptr<const IMesh> mesh,
         std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
         const FunctionSymbolId& function_id) -> std::shared_ptr<const IDiscreteFunction> {
        switch (discrete_function_descriptor->type()) {
        case DiscreteFunctionType::P0: {
          return DiscreteFunctionInterpoler{mesh, discrete_function_descriptor, function_id}.interpolate();
        }
        case DiscreteFunctionType::P0Vector: {
          return DiscreteFunctionVectorInterpoler{mesh, discrete_function_descriptor, {function_id}}.interpolate();
        }
        default: {
          throw NormalError("invalid function descriptor type");
        }
        }
      }

      ));

  this->_addBuiltinFunction("interpolate",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IDiscreteFunction>(std::shared_ptr<const IMesh>,
                                                       const std::vector<std::shared_ptr<const IZoneDescriptor>>&,
                                                       std::shared_ptr<const IDiscreteFunctionDescriptor>,
                                                       const FunctionSymbolId&)>>(
                              [](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 FunctionSymbolId& function_id) -> std::shared_ptr<const IDiscreteFunction> {
                                switch (discrete_function_descriptor->type()) {
                                case DiscreteFunctionType::P0: {
                                  return DiscreteFunctionInterpoler{mesh, interpolation_zone_list,
                                                                    discrete_function_descriptor, function_id}
                                    .interpolate();
                                }
                                case DiscreteFunctionType::P0Vector: {
                                  return DiscreteFunctionVectorInterpoler{mesh,
                                                                          interpolation_zone_list,
                                                                          discrete_function_descriptor,
                                                                          {function_id}}
                                    .interpolate();
                                }
                                default: {
                                  throw NormalError("invalid function descriptor type");
                                }
                                }
                              }

                              ));

  this
    ->_addBuiltinFunction("interpolate",
                          std::make_shared<BuiltinFunctionEmbedder<
                            std::shared_ptr<const IDiscreteFunction>(std::shared_ptr<const IMesh>,
                                                                     std::shared_ptr<const IZoneDescriptor>,
                                                                     std::shared_ptr<const IDiscreteFunctionDescriptor>,
                                                                     const FunctionSymbolId&)>>(
                            [](std::shared_ptr<const IMesh> mesh,
                               std::shared_ptr<const IZoneDescriptor> interpolation_zone,
                               std::shared_ptr<const IDiscreteFunctionDescriptor> discrete_function_descriptor,
                               const FunctionSymbolId& function_id) -> std::shared_ptr<const IDiscreteFunction> {
                              switch (discrete_function_descriptor->type()) {
                              case DiscreteFunctionType::P0: {
                                return DiscreteFunctionInterpoler{mesh,
                                                                  {interpolation_zone},
                                                                  discrete_function_descriptor,
                                                                  function_id}
                                  .interpolate();
                              }
                              case DiscreteFunctionType::P0Vector: {
                                return DiscreteFunctionVectorInterpoler{mesh,
                                                                        {interpolation_zone},
                                                                        discrete_function_descriptor,
                                                                        {function_id}}
                                  .interpolate();
                              }
                              default: {
                                throw NormalError("invalid function descriptor type");
                              }
                              }
                            }

                            ));

  this->_addBuiltinFunction("randomizeMesh",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IMesh>(std::shared_ptr<const IMesh>,
                                           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&)>>(

                              [](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::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IMesh>(std::shared_ptr<const IMesh>,
                                           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
                                           const FunctionSymbolId&)>>(

                              [](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("fixed",
                          std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryConditionDescriptor>(
                            std::shared_ptr<const IBoundaryDescriptor>)>>(

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

                            ));

  this
    ->_addBuiltinFunction("axis",
                          std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryConditionDescriptor>(
                            std::shared_ptr<const IBoundaryDescriptor>)>>(

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

                            ));

  this
    ->_addBuiltinFunction("symmetry",
                          std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryConditionDescriptor>(
                            std::shared_ptr<const IBoundaryDescriptor>)>>(

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

                            ));

  this->_addBuiltinFunction("pressure",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
                                                                  const FunctionSymbolId&)>>(

                              [](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("velocity",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
                                                                  const FunctionSymbolId&)>>(

                              [](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("dirichlet",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
                                                                  const FunctionSymbolId&)>>(

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

                              ));

  this->_addBuiltinFunction("normal_strain",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
                                                                  const FunctionSymbolId&)>>(

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

                              ));

  this->_addBuiltinFunction("fourier",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
                                                                  const FunctionSymbolId&, const FunctionSymbolId&)>>(

                              [](std::shared_ptr<const IBoundaryDescriptor> boundary, const FunctionSymbolId& alpha_id,
                                 const FunctionSymbolId& g_id) -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                return std::make_shared<FourierBoundaryConditionDescriptor>("fourier", boundary,
                                                                                            alpha_id, g_id);
                              }

                              ));

  this->_addBuiltinFunction("neumann",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
                                                                  const FunctionSymbolId&)>>(

                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
                                 const FunctionSymbolId& g_id) -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                return std::make_shared<NeumannBoundaryConditionDescriptor>("neumann", boundary, g_id);
                              }

                              ));

  this->_addBuiltinFunction("external_fsi_velocity",
                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
                                                                  const std::shared_ptr<const Socket>&)>>(

                              [](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_solver",
                            std::make_shared<BuiltinFunctionEmbedder<std::tuple<
                              std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
                              std::shared_ptr<const IDiscreteFunction>,
                              std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::vector<std::shared_ptr<
                                                                          const IBoundaryConditionDescriptor>>&,
                                                                        const double&)>>(

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

                              ));

  this->_addBuiltinFunction("eucclhyd_solver",
                            std::make_shared<BuiltinFunctionEmbedder<std::tuple<
                              std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
                              std::shared_ptr<const IDiscreteFunction>,
                              std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::vector<std::shared_ptr<
                                                                          const IBoundaryConditionDescriptor>>&,
                                                                        const double&)>>(

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

                              ));

  this->_addBuiltinFunction("heat", std::make_shared<BuiltinFunctionEmbedder<
                                      void(std::shared_ptr<const IMesh>,
                                           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
                                           const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(

                                      [](std::shared_ptr<const IMesh> p_mesh,
                                         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                           bc_descriptor_list,
                                         const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
                                         const FunctionSymbolId& f_id) -> void {
                                        switch (p_mesh->dimension()) {
                                        case 1: {
                                          HeatDiamondScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                          break;
                                        }
                                        case 2: {
                                          HeatDiamondScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                          break;
                                        }
                                        case 3: {
                                          HeatDiamondScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                          break;
                                        }
                                        default: {
                                          throw UnexpectedError("invalid mesh dimension");
                                        }
                                        }
                                      }

                                      ));

  this->_addBuiltinFunction("parabolicheat",
                            std::make_shared<BuiltinFunctionEmbedder<
                              void(std::shared_ptr<const IMesh>,
                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&,
                                   const FunctionSymbolId&, const double&, const double&)>>(

                              [](std::shared_ptr<const IMesh> p_mesh,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const FunctionSymbolId& T_id, const FunctionSymbolId& T_init_id,
                                 const FunctionSymbolId& kappa_id, const FunctionSymbolId& f_id,
                                 const double& final_time, const double& dt) -> void {
                                switch (p_mesh->dimension()) {
                                case 1: {
                                  ParabolicHeatScheme<1>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
                                                         f_id,   final_time,         dt};
                                  break;
                                }
                                case 2: {
                                  ParabolicHeatScheme<2>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
                                                         f_id,   final_time,         dt};
                                  break;
                                }
                                case 3: {
                                  ParabolicHeatScheme<3>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
                                                         f_id,   final_time,         dt};
                                  break;
                                }
                                default: {
                                  throw UnexpectedError("invalid mesh dimension");
                                }
                                }
                              }

                              ));

  this->_addBuiltinFunction(
    "parabolicheat",
    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
      const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
                               const std::shared_ptr<const IDiscreteFunction>&,
                               const std::shared_ptr<const IDiscreteFunction>&,
                               const std::shared_ptr<const IDiscreteFunction>&,
                               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&)>>(

      [](const std::shared_ptr<const IDiscreteFunction>& alpha,
         const std::shared_ptr<const IDiscreteFunction>& mub_dual,
         const std::shared_ptr<const IDiscreteFunction>& mu_dual, const std::shared_ptr<const IDiscreteFunction>& f,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
        -> const std::shared_ptr<const IDiscreteFunction> {
        return ScalarDiamondSchemeHandler{alpha, mub_dual, mu_dual, f, bc_descriptor_list}.solution();
      }

      ));
  this->_addBuiltinFunction(
    "nodalheat",
    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
      const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
                               const std::shared_ptr<const IDiscreteFunction>&,
                               const std::shared_ptr<const IDiscreteFunction>&,
                               const std::shared_ptr<const IDiscreteFunction>&,
                               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&)>>(

      [](const std::shared_ptr<const IDiscreteFunction>& alpha,
         const std::shared_ptr<const IDiscreteFunction>& mub_dual,
         const std::shared_ptr<const IDiscreteFunction>& mu_dual, const std::shared_ptr<const IDiscreteFunction>& f,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
        -> const std::shared_ptr<const IDiscreteFunction> {
        return ScalarNodalSchemeHandler{alpha, mub_dual, mu_dual, f, bc_descriptor_list}.solution();
      }

      ));

  this->_addBuiltinFunction("unsteadyelasticity",
                            std::make_shared<BuiltinFunctionEmbedder<
                              void(std::shared_ptr<const IMesh>,
                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&,
                                   const FunctionSymbolId&, const FunctionSymbolId&, const double&, const double&)>>(

                              [](std::shared_ptr<const IMesh> p_mesh,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const FunctionSymbolId& lambda_id, const FunctionSymbolId& mu_id,
                                 const FunctionSymbolId& f_id, const FunctionSymbolId& U_id,
                                 const FunctionSymbolId& U_init_id, const double& final_time,
                                 const double& dt) -> void {
                                switch (p_mesh->dimension()) {
                                case 1: {
                                  UnsteadyElasticity<1>{p_mesh, bc_descriptor_list, lambda_id,  mu_id, f_id,
                                                        U_id,   U_init_id,          final_time, dt};
                                  break;
                                }
                                case 2: {
                                  UnsteadyElasticity<2>{p_mesh, bc_descriptor_list, lambda_id,  mu_id, f_id,
                                                        U_id,   U_init_id,          final_time, dt};
                                  break;
                                }
                                case 3: {
                                  UnsteadyElasticity<3>{p_mesh, bc_descriptor_list, lambda_id,  mu_id, f_id,
                                                        U_id,   U_init_id,          final_time, dt};
                                  break;
                                }
                                default: {
                                  throw UnexpectedError("invalid mesh dimension");
                                }
                                }
                              }

                              ));

  this->_addBuiltinFunction("unsteadyelasticity",
                            std::make_shared<BuiltinFunctionEmbedder<
                              std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
                                                                       const std::shared_ptr<const IDiscreteFunction>&,
                                                                       const std::shared_ptr<const IDiscreteFunction>&,
                                                                       const std::shared_ptr<const IDiscreteFunction>&,
                                                                       const std::shared_ptr<const IDiscreteFunction>&,
                                                                       const std::shared_ptr<const IDiscreteFunction>&,
                                                                       const std::vector<std::shared_ptr<
                                                                         const IBoundaryConditionDescriptor>>&)>>(

                              [](const std::shared_ptr<const IDiscreteFunction> alpha,
                                 const std::shared_ptr<const IDiscreteFunction> lambdab,
                                 const std::shared_ptr<const IDiscreteFunction> mub,
                                 const std::shared_ptr<const IDiscreteFunction> lambda,
                                 const std::shared_ptr<const IDiscreteFunction> mu,
                                 const std::shared_ptr<const IDiscreteFunction> f,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list) -> const std::shared_ptr<const IDiscreteFunction> {
                                return VectorDiamondSchemeHandler{alpha, lambdab,           mub, lambda, mu,
                                                                  f,     bc_descriptor_list}
                                  .solution();
                              }

                              ));

  this->_addBuiltinFunction("moleculardiffusion",
                            std::make_shared<BuiltinFunctionEmbedder<std::tuple<
                              std::shared_ptr<const IDiscreteFunction>,
                              std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::vector<std::shared_ptr<
                                                                          const IBoundaryConditionDescriptor>>&)>>(
                              [](const std::shared_ptr<const IDiscreteFunction> alpha,
                                 const std::shared_ptr<const IDiscreteFunction> lambdab,
                                 const std::shared_ptr<const IDiscreteFunction> mub,
                                 const std::shared_ptr<const IDiscreteFunction> lambda,
                                 const std::shared_ptr<const IDiscreteFunction> mu,
                                 const std::shared_ptr<const IDiscreteFunction> f,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list) -> const std::tuple<std::shared_ptr<const IDiscreteFunction>,
                                                                           std::shared_ptr<const IDiscreteFunction>> {
                                return VectorDiamondSchemeHandler{alpha, lambdab,           mub, lambda, mu,
                                                                  f,     bc_descriptor_list}
                                  .apply();
                              }

                              ));

  this->_addBuiltinFunction("energybalance",
                            std::make_shared<BuiltinFunctionEmbedder<std::tuple<
                              std::shared_ptr<const IDiscreteFunction>,
                              std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::shared_ptr<const IDiscreteFunction>&,
                                                                        const std::vector<std::shared_ptr<
                                                                          const IBoundaryConditionDescriptor>>&)>>(
                              [](const std::shared_ptr<const IDiscreteFunction> lambdab,
                                 const std::shared_ptr<const IDiscreteFunction> mub,
                                 const std::shared_ptr<const IDiscreteFunction> U,
                                 const std::shared_ptr<const IDiscreteFunction> dual_U,
                                 const std::shared_ptr<const IDiscreteFunction> source,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list) -> const std::tuple<std::shared_ptr<const IDiscreteFunction>,
                                                                           std::shared_ptr<const IDiscreteFunction>> {
                                return EnergyComputerHandler{lambdab, mub, U, dual_U, source, bc_descriptor_list}
                                  .computeEnergyUpdate();
                              }

                              ));

  this->_addBuiltinFunction("heat2",
                            std::make_shared<BuiltinFunctionEmbedder<
                              void(std::shared_ptr<const IMesh>,
                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(

                              [](std::shared_ptr<const IMesh> p_mesh,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
                                 const FunctionSymbolId& f_id) -> void {
                                switch (p_mesh->dimension()) {
                                case 1: {
                                  HeatDiamondScheme2<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                  break;
                                }
                                case 2: {
                                  HeatDiamondScheme2<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                  break;
                                }
                                case 3: {
                                  HeatDiamondScheme2<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                  break;
                                }
                                default: {
                                  throw UnexpectedError("invalid mesh dimension");
                                }
                                }
                              }

                              ));

  this->_addBuiltinFunction("elasticity",
                            std::make_shared<BuiltinFunctionEmbedder<
                              void(std::shared_ptr<const IMesh>,
                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&,
                                   const FunctionSymbolId&)>>(

                              [](std::shared_ptr<const IMesh> p_mesh,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const FunctionSymbolId& lambda_id, const FunctionSymbolId& mu_id,
                                 const FunctionSymbolId& f_id, const FunctionSymbolId& U_id) -> void {
                                switch (p_mesh->dimension()) {
                                case 1: {
                                  ElasticityDiamondScheme<1>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
                                  break;
                                }
                                case 2: {
                                  ElasticityDiamondScheme<2>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
                                  break;
                                }
                                case 3: {
                                  ElasticityDiamondScheme<3>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
                                  break;
                                }
                                default: {
                                  throw UnexpectedError("invalid mesh dimension");
                                }
                                }
                              }

                              ));

  this->_addBuiltinFunction("heat5points",
                            std::make_shared<BuiltinFunctionEmbedder<
                              void(std::shared_ptr<const IMesh>,
                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId)>>(

                              [](std::shared_ptr<const IMesh> p_mesh,
                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                   bc_descriptor_list,
                                 const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
                                 const FunctionSymbolId& f_id) -> void {
                                switch (p_mesh->dimension()) {
                                case 1: {
                                  Heat5PointsAlgorithm<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                  break;
                                }
                                case 2: {
                                  Heat5PointsAlgorithm<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                  break;
                                }
                                case 3: {
                                  Heat5PointsAlgorithm<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                  break;
                                }
                                default: {
                                  throw UnexpectedError("invalid mesh dimension");
                                }
                                }
                              }

                              ));

  this
    ->_addBuiltinFunction("lagrangian",
                          std::make_shared<BuiltinFunctionEmbedder<
                            std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IMesh>&,
                                                                     const std::shared_ptr<const IDiscreteFunction>&)>>(

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

                            ));

  this->_addBuiltinFunction("acoustic_dt",
                            std::make_shared<
                              BuiltinFunctionEmbedder<double(const std::shared_ptr<const IDiscreteFunction>&)>>(

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

                              ));

  this
    ->_addBuiltinFunction("cell_volume",
                          std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IDiscreteFunction>(
                            const std::shared_ptr<const IMesh>&)>>(

                            [](const std::shared_ptr<const IMesh>& i_mesh) -> std::shared_ptr<const IDiscreteFunction> {
                              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<const DiscreteFunctionP0<
                                  Dimension, double>>(mesh, copy(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<const DiscreteFunctionP0<
                                  Dimension, double>>(mesh, copy(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<const DiscreteFunctionP0<
                                  Dimension, double>>(mesh, copy(MeshDataManager::instance().getMeshData(*mesh).Vj()));
                              }
                              default: {
                                throw UnexpectedError("invalid mesh dimension");
                              }
                              }
                            }

                            ));

  MathFunctionRegisterForVh{*this};
}

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