#ifndef HYPERELASTIC_SOLVER_HPP
#define HYPERELASTIC_SOLVER_HPP

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

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

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

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

 private:
  struct IHyperelasticSolver
  {
    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::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>& CG,
                 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>& CG,
          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;

    IHyperelasticSolver()                                 = default;
    IHyperelasticSolver(IHyperelasticSolver&&)            = default;
    IHyperelasticSolver& operator=(IHyperelasticSolver&&) = default;

    virtual ~IHyperelasticSolver() = default;
  };

  template <size_t Dimension>
  class HyperelasticSolver;

  std::unique_ptr<IHyperelasticSolver> m_hyperelastic_solver;

 public:
  const IHyperelasticSolver&
  solver() const
  {
    return *m_hyperelastic_solver;
  }

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

#endif   // HYPERELASTIC_SOLVER_HPP