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

#include <dolfin.h>

#include "SubdomainNumbering.h"
#include "RFASolver.h"

using namespace dolfin;

constexpr real radius = 0.001155;
constexpr real xc = 0.04;
constexpr real yc = 0.04;

constexpr real density_b = 1050;
constexpr real mu = 2.65e-3;
constexpr real nu = mu/density_b;

constexpr real bodyTemp = 37 + 273.15;

constexpr real cfl_target = 100;

// Inflow velocity
class SalineInflow : public Function
{
public:
    SalineInflow(Mesh& mesh, real u) : Function(mesh), u(u) {}

    void eval(real* values, const real* x) const
    {
        values[0] = u*(x[0]-xc)/radius;
        values[1] = u*(x[1]-yc)/radius;
        values[2] = 0.0;
    }

private:
    real u;
};

int main(int argc, char* argv[])
{
    dolfin_init(argc, argv);
    {
        if(!ParameterSystem::parameters.defined("adapt_percent"))
            dolfin_add("adapt_percent", 10.0);
        if(!ParameterSystem::parameters.defined("alpha"))
            dolfin_add("alpha", 0.0);

        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 ablationTime;
        real nsInitTime;
        real Uin;
        real density_t;
        real specificHeat0_t;
        real specificHeat1_t;
        real kappa0_t;
        real kappa1_t;
        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")
            ("saline_inflow_modulus",
                 po::value<real>(&salineInflowModulus),
                 "Modulus of the saline inflow velocity")
            ("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")
            ("ablationTime",
                 po::value<real>(&ablationTime),
                 "Final time for the ablation")
            ("ns_init_time",
                po::value<real>(&nsInitTime),
                "Initialization time for Navier-Stokes")
            ("ns_inflow_modulus",
                 po::value<real>(&Uin),
                 "Inflow Navier-Stokes")
            ("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")
            ("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")
            ("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).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,
            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);
        //  dolfin_set("Krylov monitor convergence",true);

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

        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 catheterBodyMesh((mesh_prefix + "CatheterBody" + 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 inletMesh((mesh_prefix + "Inlet" + mesh_suffix).c_str());
        Mesh interfaceMesh((mesh_prefix + "Interface" + 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 thermistorUpperInterfaceMesh((mesh_prefix + "ThermistorUpperInterface" + mesh_suffix).c_str());
        Mesh topMesh((mesh_prefix + "Top" + mesh_suffix).c_str());
        Mesh wallsDownstreamMesh((mesh_prefix + "WallsDownstream" + mesh_suffix).c_str());

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

        SubDomain boardDomain(boardDomainMesh);
        SubDomain bottom(bottomMesh);
        SubDomain catheterBody(catheterBodyMesh);
        SubDomain fluidDomain(fluidDomainMesh);
        SubDomain electrode(electrodeMesh);
        SubDomain electrodeUpperInterface(electrodeUpperInterfaceMesh);
        SubDomain inlet(inletMesh);
        SubDomain interface(interfaceMesh);
        SubDomain leftFluidWall(leftFluidWallMesh);
        SubDomain outlet(outletMesh);
        SubDomain rightFluidWall(rightFluidWallMesh);
        SubDomain salineInlet(salineInletMesh);
        SubDomain solidDomain(solidDomainMesh);
        SubDomain solidWalls(solidWallsMesh);
        SubDomain thermistor(thermistorMesh);
        SubDomain thermistorUpperInterface(thermistorUpperInterfaceMesh);
        SubDomain top(topMesh);
        SubDomain wallsDownstream(wallsDownstreamMesh);

        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 Navier-Stokes
        Function inflow(mesh,Array<real>(Uin,0,0));
        SalineInflow salineInflow(mesh,salineInflowModulus);
        Function noSlipVelocity(mesh, Array<real>(0,0,0));
        Function outflow(mesh,0);
        Function u_i0(mesh, Array<real>(0,0,0));
        Function f(mesh, Array<real>(0,0,0));

        MeshFunction<uint> thetadrag(mesh,mesh.topology().dim()-1);
        thetadrag = 1;
        Lift lift(mesh,&thetadrag);
        OutputQuantityList outputQuantities;
        outputQuantities.push_back(&lift);
        MeshFunction<uint> sm(mesh,mesh.topology().dim()-1);
        sm = 1;

        DirichletBCList dbcs_m;
        dbcs_m.push_back(std::make_pair(&inlet,&inflow));
        dbcs_m.push_back(std::make_pair(&salineInlet,&salineInflow));
        dbcs_m.push_back(std::make_pair(&catheterBody,&noSlipVelocity));
        dbcs_m.push_back(std::make_pair(&top,&noSlipVelocity));
        dbcs_m.push_back(std::make_pair(&interface,&noSlipVelocity));
        dbcs_m.push_back(std::make_pair(&leftFluidWall,&noSlipVelocity));
        dbcs_m.push_back(std::make_pair(&rightFluidWall,&noSlipVelocity));
        dbcs_m.push_back(std::make_pair(&solidWalls,&noSlipVelocity));
        dbcs_m.push_back(std::make_pair(&bottom,&noSlipVelocity));
        dbcs_m.push_back(std::make_pair(&solidDomain,&noSlipVelocity));
        dbcs_m.push_back(std::make_pair(&boardDomain,&noSlipVelocity));
        dbcs_m.push_back(std::make_pair(&electrode,&noSlipVelocity));
        dbcs_m.push_back(std::make_pair(&thermistor,&noSlipVelocity));

        DirichletBCList dbcs_c;
        dbcs_c.push_back(std::make_pair(&wallsDownstream,&outflow));
        //  dbcs_c.push_back(std::make_pair(&solidDomain,&outflow));

        DirichletBCList dbcs_dm;

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

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

        // Parameters for Bioheat
        Function boundaryTemperature(mesh,310.15);
        Function tInitial(mesh,bodyTemp);
        Function ts(mesh,0);
        Function salineTemperature(mesh,295.15);

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

        RFASolver rfaSolver(
                mesh,
                subdomains,
                dbcs_m,
                dbcs_c,
                dbcs_dm,
                NULL,
                &u_i0,
                outputQuantities,
                &f,
                cfl_target,
                nu,
                &ps,
                &electrodeUpperInterface,
                dbcs_p,
                &ts,
                &tInitial,
                dbcs_b,
                tissueParameters);
        rfaSolver.setDt(dt);
        rfaSolver.setNoSamples(no_samples);
        rfaSolver.setTFinal(tFinal);
        rfaSolver.setAblationTime(ablationTime);
        rfaSolver.setNSInitTime(nsInitTime);
        rfaSolver.setTargetPower(targetPower);
        rfaSolver.setTissuePowerFraction(tissuePowerFraction);
        rfaSolver.setInitialResistance(initialResistance);
        rfaSolver.setSigmaBoard(0.123);
        rfaSolver.setMinSigmaBoard(0);
        rfaSolver.setMaxSigmaBoard(rfaSolver.getSigmaBoard()*2);
        rfaSolver.setPowerTolerance(1e-2);

        rfaSolver.run();
    }

    dolfin_finalize();
    return 0;
}
