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

#include <dolfin.h>

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

using namespace dolfin;

int main(int argc, char* argv[])
{
{
    std::string config_file;
    std::string mesh_prefix;
    std::string mesh_suffix;
    real targetPower;
    real tissuePowerFraction;
    real initialResistance;
    real salineInflowModulus;
    real dt;
    int no_samples;
    real tFinal;
    real nsInitTime;
    real sigma0_t;
    real sigma1_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()
        ("targetPower",
         po::value<real>(&targetPower),
         "Catheter targetPower")
        ("tissuePowerFraction",
         po::value<real>(&tissuePowerFraction),
         "Fraction of the power supplied to the tissue")
        ("initialResistance",
         po::value<real>(&initialResistance),
         "Initial resistance")
        ("mesh_prefix",
         po::value<std::string>(&mesh_prefix)->default_value(std::string()))
        ("mesh_suffix",
         po::value<std::string>(&mesh_suffix)->default_value(std::string()))
        ("sigma0_t",
         po::value<real>(&sigma0_t),
         "Sigma of the tissue, first parameter")
        ("sigma1_t",
         po::value<real>(&sigma1_t),
         "Sigma 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).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(
            sigma0_t,
            sigma1_t,
            0,
            0,
            0,
            0,
            0);

    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());
    //  Mesh mesh((mesh_prefix + "Domain_elastic_big_10g.bin").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 electrodeUpperInterfaceMesh((mesh_prefix + "ElectrodeUpperInterface" + mesh_suffix).c_str());
    Mesh solidDomainMesh((mesh_prefix + "SolidDomain" + mesh_suffix).c_str());
    Mesh thermistorMesh((mesh_prefix + "Thermistor" + mesh_suffix).c_str());

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

    SubDomain boardDomain(boardDomainMesh);
    SubDomain bottom(bottomMesh);
    SubDomain fluidDomain(fluidDomainMesh);
    SubDomain electrode(electrodeMesh);
    SubDomain electrodeUpperInterface(electrodeUpperInterfaceMesh);
    SubDomain solidDomain(solidDomainMesh);
    SubDomain thermistor(thermistorMesh);

    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);

    // Parameters for potential
    Function boundaryPotential(mesh,0);
    Function potentialSource(mesh,0);

    DirichletBCList dbcs_p;
    dbcs_p.push_back(std::make_pair(&bottom,&boundaryPotential));

    PotentialProblem potentialProblem(
            mesh,
            subdomains,
            &potentialSource,
            &electrode,
            dbcs_p,
            tissueParameters);
    potentialProblem.setTargetPower(targetPower);
    potentialProblem.setTissuePowerFraction(tissuePowerFraction);
    potentialProblem.setInitialResistance(initialResistance);

    PotentialSolver potentialSolver(&potentialProblem);
    Function bodyTemperature(mesh,273.15+37);
    potentialSolver.init(&bodyTemperature);
    potentialSolver.setSigmaBoard(0.123);
    potentialSolver.setMinSigmaBoard(0);
    potentialSolver.setMaxSigmaBoard(potentialSolver.getSigmaBoard()*2);
    potentialSolver.setPowerTolerance(1e-8);

    potentialSolver.solve();
    potentialSolver.save();
}

    dolfin_finalize();
    return 0;
}
