#include <boost/program_options.hpp>
namespace po = boost::program_options;

#include <dolfin.h>

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

using namespace dolfin;

const real density = 1;
const real specificHeat = 1;
const real bodyTemp = 37 + 273.15;

class Boundary : public SubDomain
{
    bool inside(const real* x, bool on_boundary) const
    {
        return on_boundary;
    }
};

int main(int argc, char* argv[])
{
    dolfin_init(argc,argv);

    std::string config_file;
    std::string mesh_prefix;
    std::string mesh_suffix;
    real dt;
    int no_samples;
    real tFinal;
    real density_t;
    real specificHeat0_t;
    real specificHeat1_t;
    real kappa0_t;
    real kappa1_t;

    po::options_description generic("Generic options");
    generic.add_options()
        ("config,c", po::value<std::string>(&config_file),
         "Configuration file");

    po::options_description config("Configuration");
    config.add_options()
        ("no_samples",
         po::value<int>(&no_samples),
         "Number of samples to output")
        ("dt",
         po::value<real>(&dt),
         "Time step for bioheat equation")
        ("tFinal",
         po::value<real>(&tFinal),
         "Final time for the simulation")
        ("mesh_prefix",
         po::value<std::string>(&mesh_prefix)->default_value(std::string()))
        ("mesh_suffix",
         po::value<std::string>(&mesh_suffix)->default_value(std::string()))
        ("density_t",
         po::value<real>(&density_t),
         "Density of the tissue")
        ("kappa0_t",
         po::value<real>(&kappa0_t),
         "Kappa of the tissue, first parameter")
        ("kappa1_t",
         po::value<real>(&kappa1_t),
         "Kappa of the tissue, second parameter")
        ("specificHeat0_t",
         po::value<real>(&specificHeat0_t),
         "Specific heat of the tissue, first parameter")
        ("specificHeat1_t",
         po::value<real>(&specificHeat1_t),
         "Specific heat of the tissue, second parameter")
        ;

    po::options_description cmdline_options;
    cmdline_options.add(generic);

    po::options_description config_file_options;
    config_file_options.add(config);

    po::variables_map vm;
    store(po::command_line_parser(argc, argv).
            options(cmdline_options).allow_unregistered().run(), vm);
    notify(vm);

    std::ifstream ifs(config_file.c_str());
    if (!ifs)
    {
        std::cout << "can not open config file: " << config_file << "\n";
        return 0;
    }
    else
    {
        store(parse_config_file(ifs, config_file_options), vm);
        notify(vm);
    }

    TissueParameters tissueParameters(
            0,
            0,
            kappa0_t,
            kappa1_t,
            specificHeat0_t,
            specificHeat1_t,
            density_t);

    dolfin_set("output destination","silent");
    if(dolfin::MPI::processNumber() == 0)
        dolfin_set("output destination","terminal");

    dolfin_set("SubDomain Intersect Boundary",false);

    // Create mesh
    Mesh mesh((mesh_prefix + "Domain" + mesh_suffix).c_str());

    File meshfile("mesh_out.bin");
    meshfile << mesh;

    dolfin_set("Mesh read in serial", true);

    // Create boundary conditions
    Mesh boardDomainMesh((mesh_prefix + "BoardDomain" + mesh_suffix).c_str());
    Mesh bottomMesh((mesh_prefix + "Bottom" + mesh_suffix).c_str());
    Mesh fluidDomainMesh((mesh_prefix + "FluidDomain" + mesh_suffix).c_str());
    Mesh electrodeMesh((mesh_prefix + "Electrode" + mesh_suffix).c_str());
    Mesh inletMesh((mesh_prefix + "Inlet" + mesh_suffix).c_str());
    Mesh leftFluidWallMesh((mesh_prefix + "LeftFluidWall" + mesh_suffix).c_str());
    Mesh outletMesh((mesh_prefix + "Outlet" + mesh_suffix).c_str());
    Mesh rightFluidWallMesh((mesh_prefix + "RightFluidWall" + mesh_suffix).c_str());
    Mesh salineInletMesh((mesh_prefix + "SalineInlet" + mesh_suffix).c_str());
    Mesh solidDomainMesh((mesh_prefix + "SolidDomain" + mesh_suffix).c_str());
    Mesh solidWallsMesh((mesh_prefix + "SolidWalls" + mesh_suffix).c_str());
    Mesh thermistorMesh((mesh_prefix + "Thermistor" + mesh_suffix).c_str());
    Mesh topMesh((mesh_prefix + "Top" + mesh_suffix).c_str());

    dolfin_set("Mesh read in serial", false);

    SubDomain boardDomain(boardDomainMesh);
    SubDomain bottom(bottomMesh);
    SubDomain fluidDomain(fluidDomainMesh);
    SubDomain electrode(electrodeMesh);
    SubDomain inlet(inletMesh);
    SubDomain leftFluidWall(leftFluidWallMesh);
    SubDomain outlet(outletMesh);
    SubDomain rightFluidWall(rightFluidWallMesh);
    SubDomain salineInlet(salineInletMesh);
    SubDomain solidDomain(solidDomainMesh);
    SubDomain solidWalls(solidWallsMesh);
    SubDomain thermistor(thermistorMesh);
    SubDomain top(topMesh);

    std::vector<SubDomain*> subdomains(5);
    subdomains[i_t] = (&solidDomain);
    subdomains[i_b] = (&fluidDomain);
    subdomains[i_e] = (&electrode);
    subdomains[i_th] = (&thermistor);
    subdomains[i_bo] = (&boardDomain);
    Function catheterTemperature(mesh,1);

    Boundary boundary;

    Function psi(mesh,0);
    Function vel(mesh,Array<real>(0,0,0));
    Function boundaryTemperature(mesh,310.15);
    Function tInitial(mesh,bodyTemp);
    Function salineTemperature(mesh,295.15);

    DirichletBCList dbcs;
    dbcs.push_back(std::make_pair(&leftFluidWall,&boundaryTemperature));
    dbcs.push_back(std::make_pair(&rightFluidWall,&boundaryTemperature));
    dbcs.push_back(std::make_pair(&solidWalls,&boundaryTemperature));
    dbcs.push_back(std::make_pair(&top,&boundaryTemperature));
    dbcs.push_back(std::make_pair(&bottom,&boundaryTemperature));
    dbcs.push_back(std::make_pair(&inlet,&boundaryTemperature));
    dbcs.push_back(std::make_pair(&outlet,&boundaryTemperature));
    dbcs.push_back(std::make_pair(&salineInlet,&salineTemperature));

    BioProblem problem(
            mesh,
            subdomains,
            &vel,
            &psi,
            &tInitial,
            dbcs,
            tissueParameters);
    BioSolver solver(&problem);
    solver.setTFinal(1e-1);

    solver.run();

    dolfin_finalize();

    return 0;
}
