#include "BioSolver.h"
#include "SubdomainNumbering.h"

BioSolver::BioSolver(BioProblem* problem) :
    problem(problem),
    no_samples(100),
    sample(1)
{
    Mesh& mesh = problem->getMesh();

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

#warning "This should be fixed"
    T = new Function(mesh,37+273.15);
    TOld = new Function(mesh,37+273.15);
    //  *T = *problem->getTInitial();
    //  *TOld = *problem->getTInitial();
    dt = 0.01;
    dt_f = new Function(mesh,dt);
    h = new MeshSize(mesh);
    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);

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

    a = new BioheatBilinearForm(mesh,forms);
    (*a)["vel"] = problem->getVelocity();
    (*a)["density_t"] = problem->getTissueDensity();
    (*a)["specificHeat0_t"] = problem->getTissueSpecificHeat0();
    (*a)["specificHeat1_t"] = problem->getTissueSpecificHeat1();
    (*a)["kappa0_t"] = problem->getTissueKappa0();
    (*a)["kappa1_t"] = problem->getTissueKappa1();
    (*a)["uOld"] = TOld;
    (*a)["dt"] = dt_f;
    (*a)["h"] = h;

    L = new BioheatLinearForm(mesh,forms);
    (*L)["Phi"] = problem->getPhi();
    (*L)["sigma_bo"] = sigmaBoard;
    (*L)["density_t"] = problem->getTissueDensity();
    (*L)["specificHeat0_t"] = problem->getTissueSpecificHeat0();
    (*L)["specificHeat1_t"] = problem->getTissueSpecificHeat1();
    (*L)["sigma0_t"] = problem->getTissueSigma0();
    (*L)["sigma1_t"] = problem->getTissueSigma1();

    pde = new LinearPDE(*a, *L, mesh, bcs);
    pde->setCellDomains(domains);

    tInitial = 0;
    t = tInitial;
    tFinal = 1;

    file_T = new File("temperature.bin");
}

void BioSolver::solveSystem()
{
    pde->solve(*T);
}

real BioSolver::step()
{
    solveSystem();

    t += dt;

    if(t >= tFinal*(real(sample)/real(no_samples)))
    {
        cout << "saving bioheat at time " << t << endl;
        save();
        sample++;
    }

    *TOld = *T;

    return t;
}

void BioSolver::run()
{
    while(t <= tFinal)
    {
        step();
    }
}

void BioSolver::save()
{
    *file_T << *T;
}

BioSolver::~BioSolver ()
{
    bcs.clear();
    delete T;
    delete TOld;
    delete a;
    delete L;
    delete pde;
    delete file_T;
    delete h;
    delete dt_f;
    delete domains;
    delete sigmaBoard;
}
