#ifndef HYPOELASTIC_SOLVER_HPP
#define HYPOELASTIC_SOLVER_HPP

#include <memory>
#include <tuple>
#include <vector>

class IBoundaryConditionDescriptor;
class IMesh;
class ItemValueVariant;
class SubItemValuePerItemVariant;
class DiscreteFunctionVariant;

double hypoelastic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c);

class HypoelasticSolverHandler
{
 public:
  enum class SolverType
  {
    Glace,
    Eucclhyd
  };

 private:
  struct IHypoelasticSolver
  {
    virtual std::tuple<const std::shared_ptr<const ItemValueVariant>,
                       const std::shared_ptr<const SubItemValuePerItemVariant>>
    compute_fluxes(
      const SolverType& solver_type,
      const std::shared_ptr<const DiscreteFunctionVariant>& rho,
      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::shared_ptr<const DiscreteFunctionVariant>& p,
      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;

    virtual 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>>
    apply_fluxes(const double& dt,
                 const std::shared_ptr<const DiscreteFunctionVariant>& rho,
                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
                 const std::shared_ptr<const ItemValueVariant>& ur,
                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const = 0;

    virtual 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>>
    apply(const SolverType& solver_type,
          const double& dt,
          const std::shared_ptr<const DiscreteFunctionVariant>& rho,
          const std::shared_ptr<const DiscreteFunctionVariant>& u,
          const std::shared_ptr<const DiscreteFunctionVariant>& E,
          const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
          const std::shared_ptr<const DiscreteFunctionVariant>& aL,
          const std::shared_ptr<const DiscreteFunctionVariant>& aT,
          const std::shared_ptr<const DiscreteFunctionVariant>& p,
          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;

    IHypoelasticSolver()                                = default;
    IHypoelasticSolver(IHypoelasticSolver&&)            = default;
    IHypoelasticSolver& operator=(IHypoelasticSolver&&) = default;

    virtual ~IHypoelasticSolver() = default;
  };

  template <size_t Dimension>
  class HypoelasticSolver;

  std::unique_ptr<IHypoelasticSolver> m_hypoelastic_solver;

 public:
  const IHypoelasticSolver&
  solver() const
  {
    return *m_hypoelastic_solver;
  }

  HypoelasticSolverHandler(const std::shared_ptr<const IMesh>& mesh);
};

#endif   // HYPOELASTIC_SOLVER_HPP