#include "PotentialSolver.h"
#include "SubdomainNumbering.h"

PotentialSolver::PotentialSolver(PotentialProblem* problem) :
    problem(problem),
    powerTolerance(0.1),
    a(0),
    L(0),
    pde(0),
    M(0),
    minSigmaBoard(0),
    maxSigmaBoard(10),
    domainPower(15)
{
    Mesh& mesh = problem->getMesh();

    domains = new MeshFunction<uint>(mesh,mesh.topology().dim());
    *domains = 99;
    const std::vector<SubDomain*>& sd = problem->getSubdomains();
    sd[i_e]->mark(*domains,i_e);
    sd[i_th]->mark(*domains,i_th);
    sd[i_bo]->mark(*domains,i_bo);
    sd[i_t]->mark(*domains,i_t);
    sd[i_b]->mark(*domains,i_b);

    File subfile("subdomains.bin");
    subfile << *domains;

    Phi = new Function;
    assembler = new Assembler(mesh);

    electrodePotential = new ElectrodePotential(problem->getMesh());
    sigmaBoard = new SigmaBoard(problem->getMesh());

    bcs.push_back(new DirichletBC(*electrodePotential,
                  problem->getMesh(),
                  *problem->getElectrode()));
    const DirichletBCList& dirbcs = problem->getDirichletBCs();
    for (int i = 0; i < dirbcs.size(); i++)
    {
        bcs.push_back(new DirichletBC(*dirbcs[i].second,
                      problem->getMesh(),
                      *dirbcs[i].first));
    }

    setPotential(1);
}

#warning "Why not just assigning?"
void PotentialSolver::init(Function* T)
{
    a = new PotentialBilinearForm(problem->getMesh(),forms);
    (*a)["T"] = T;
    (*a)["sigma_bo"] = sigmaBoard;
    (*a)["sigma0_t"] = problem->getTissueSigma0();
    (*a)["sigma1_t"] = problem->getTissueSigma1();

    L = new PotentialLinearForm(problem->getMesh(),forms);
    (*L)["f"] = problem->getSource();
    Phi->init(problem->getMesh(),*a,0);
    M = new PotentialFunctional(problem->getMesh(),forms);
    (*M)["unum"] = Phi;
    (*M)["sigma0_t"] = problem->getTissueSigma0();
    (*M)["sigma1_t"] = problem->getTissueSigma1();
    pde = new LinearPDE(*a, *L, problem->getMesh(), bcs);
    pde->setCellDomains(domains);

    file_Phi = new File("potential.bin");
}

PotentialSolver::~PotentialSolver()
{
    bcs.clear();
    delete Phi;
    delete a;
    delete L;
    delete M;
    delete assembler;
    delete pde;
    delete file_Phi;
    delete domains;
    delete electrodePotential;
    delete sigmaBoard;
}

void PotentialSolver::solve()
{
    real powerDelta = getPowerTolerance() + 1;

    while (powerDelta > getPowerTolerance())
    {
        pde->solve(*Phi);

        real newPower = assembler->assemble(
                *M,
                *domains,
                MeshFunction<uint>(),
                MeshFunction<uint>());
        cout << "Current power: " << newPower << endl;
        cout << "Target power: " << domainPower << endl;
        powerDelta = fabs(domainPower-newPower);

        setPotential(getPotential()*sqrt(domainPower/newPower));
    }
}

void PotentialSolver::save()
{
    *file_Phi << *Phi;
}

void PotentialSolver::matchInitialPower()
{
    cout << "Matching tissue power of " << problem->getTargetPower()*problem->getTissuePowerFraction() << "W" << endl;
    setPotential(sqrt(problem->getInitialResistance()*problem->getTargetPower()));
    cout << "Potential set to  " << getPotential() << "V" << endl;

    real powerDelta = getPowerTolerance() + 1;

    Mesh& mesh = problem->getMesh();
    MeshFunction<uint> tissuePowerDomains(mesh,mesh.topology().dim());
    const std::vector<SubDomain*>& sd = problem->getSubdomains();

    tissuePowerDomains = 99;
    sd[i_t]->mark(tissuePowerDomains,i_t);

    while (powerDelta > getPowerTolerance())
    {
        setSigmaBoard(0.5*(getMinSigmaBoard()+getMaxSigmaBoard()));
        pde->solve(*Phi);

        domainPower = assembler->assemble(
                *M,
                *domains,
                MeshFunction<uint>(),
                MeshFunction<uint>());
        cout << "Current sigma: " << getSigmaBoard() << endl;
        cout << "Current power: " << domainPower << endl;

        real tissuePower = assembler->assemble(
                *M,
                tissuePowerDomains,
                MeshFunction<uint>(),
                MeshFunction<uint>());
        cout << "Current tissue power: " << tissuePower << endl;

        powerDelta = fabs(problem->getTargetPower()*problem->getTissuePowerFraction()-tissuePower);

        if (tissuePower > problem->getTargetPower()*problem->getTissuePowerFraction())
        {
            setMaxSigmaBoard(getSigmaBoard());
        }
        else
        {
            setMinSigmaBoard(getSigmaBoard());
        }
    }
}
