#ifndef POLYSOLVER_H
#define POLYSOLVER_H

#include "mfem.hpp"
#include <iostream>
#include <string>
#include <memory>

#include "meshIO.h"
#include "polyCoeff.h"
#include "polyMFEMUtils.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 thetaSerieseExpansion(double xi, double n, int order);
}

class PolySolver {
private:
    Config& config = Config::getInstance();
    Probe::LogManager& logManager = Probe::LogManager::getInstance();
    quill::Logger* logger;
    double n, order;
    MeshIO meshIO;
    mfem::Mesh& mesh;
    
    std::unique_ptr<mfem::H1_FECollection> feCollection;

    std::unique_ptr<mfem::FiniteElementSpace> feSpace;

    std::unique_ptr<polyMFEMUtils::CompositeNonlinearIntegrator> compositeIntegrator;

    std::unique_ptr<mfem::NonlinearForm> nonlinearForm;
    std::unique_ptr<mfem::LinearForm> C; // For the constraint equation

    std::unique_ptr<mfem::GridFunction> u;

    std::unique_ptr<mfem::VectorConstantCoefficient> diffusionCoeff;
    std::unique_ptr<mfem::ConstantCoefficient> nonLinearSourceCoeff;
    std::unique_ptr<polyMFEMUtils::GaussianCoefficient> gaussianCoeff;

    double C_val;

    void assembleNonlinearForm();
    void assembleConstraintForm();


public:
    PolySolver(double n, double order);
    ~PolySolver();
    void solve();
    mfem::Mesh& getMesh() { return mesh; }
    mfem::GridFunction& getSolution() { return *u; }
    double getN() { return n; }
    double getOrder() { return order; }
};

#endif // POLYSOLVER_H