#pragma once

#include "mfem.hpp"
#include <memory>
#include <utility>

#include "integrators.h"
#include "4DSTARTypes.h"
#include "operator.h"
#include "config.h"
#include "probe.h"
#include "quill/Logger.h"


namespace laneEmden {
    double a (int k, double n);
    double c(int m, double n);
    double thetaSeriesExpansion(double xi, double n, int order);
}

// Struct to persist lifetime of the linear and nonlinear solvers
struct solverBundle {
    mfem::GMRESSolver solver; // Must be first so it lives longer than the newton solver
    mfem::NewtonSolver newton; // Must be second so that when it is destroyed the solver is still alive preventing a double delete
};

class PolySolver {
public: // Public methods
    PolySolver(double n, double order);
    ~PolySolver();

    void solve() const;

    double getN() const { return m_polytropicIndex; }
    double getOrder() const { return m_feOrder; }
    mfem::Mesh* getMesh() const { return m_mesh.get(); }
    mfem::GridFunction& getSolution() const { return *m_theta; }

private: // Private Attributes
    Config& m_config = Config::getInstance();
    Probe::LogManager& m_logManager = Probe::LogManager::getInstance();
    quill::Logger* m_logger = m_logManager.getLogger("log");
    double m_polytropicIndex, m_feOrder;
    
    std::unique_ptr<mfem::Mesh> m_mesh;
    std::unique_ptr<mfem::H1_FECollection> m_fecH1;
    std::unique_ptr<mfem::RT_FECollection> m_fecRT;

    std::unique_ptr<mfem::FiniteElementSpace> m_feTheta;
    std::unique_ptr<mfem::FiniteElementSpace> m_fePhi;

    std::unique_ptr<mfem::GridFunction> m_theta;
    std::unique_ptr<mfem::GridFunction> m_phi;

    std::unique_ptr<PolytropeOperator> m_polytropOperator;

    std::unique_ptr<mfem::OperatorJacobiSmoother> m_prec;


private: // Private methods
    void assembleBlockSystem();
    SSE::MFEMArrayPairSet getEssentialTrueDof() const;
    std::pair<mfem::Array<int>, mfem::Array<int>> findCenterElement() const;
    void setInitialGuess() const;
    void saveAndViewSolution(const mfem::BlockVector& state_vector) const;
    solverBundle setupNewtonSolver() const;
    void setupOperator() const;

    void LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel,
                              double &gmresRelTol, double &gmresAbsTol, int &gmresMaxIter, int &gmresPrintLevel) const;

};