#include "geom.h"
#include "getmatrix.h"
#include "jacsolve.h"
#include "userset.h"
#include "vdop.h"

#include <cassert>
#include <chrono>           // timing
#include <cmath>
#include <iostream>
#include <omp.h>

#include "precice/precice.hpp"
#include <sstream>
using namespace precice;

using namespace std;
using namespace std::chrono;  // timing


int main(int argc, char **argv )
{
    // generating the mesh
    Mesh const mesh_c("../generate_mesh/coffee_cup.txt", "../generate_mesh/coffee_cup_sd.txt");
    bool ba = mesh_c.checkObtuseAngles();
    if (ba) cout << "mesh corrected" << endl;
    //mesh_c.DebugEdgeBased();

    //gMesh_Hierarchy ggm(mesh_c, nrefine);
    //const Mesh &mesh = ggm.finest();
    //mesh.Debug();
    //mesh_c.DebugEdgeBased();

    // ##########################################
    // Parameteres
    // ##########################################

    double dt = 1.0;           // time step
    int steps = 100;            // number of time iterations

    double u0_mug = 18.0;
    double u0_fluid = 80.0;
    double u0_air = 18.0;
    double u_out = 18.0;

    // ##########################################
    // Assembling
    // ##########################################

    // Initializing FEM matrix !pattern! (only zero entries now)
    FEM_Matrix SK(mesh_c);                                    // CRS matrix
    //SK.writeBinary("sparseMatrix.bin");
    //SK.Debug();

    vector<double> fv(SK.Nrows(), 0.0);                    
    SK.CalculateRHS(fv, [](double x, double y) {return 0;});  // rhs (f=0)

    SK.CalculateLaplace_mult(fv);                             // stiffness matrix (+K)
    SK.AddMass_mult(fv, 1.0/dt);                              // add mass matrix (+M/dt)
    SK.ApplyRobinBC_mult(fv, u_out);                           // apply Robin-bnd (+C = +F)
                                                              // SK = M/dt + K + C = F
    
    // SK.Debug();
    // SK.CheckRowSum();
    // SK.CheckMatrix();
    
    FEM_Matrix Mdt(mesh_c);
    Mdt.AddMass_mult(fv, 1.0/dt);                              // mass matrix (Mdt = M/dt)

    // Mdt.Debug();
    // Mdt.CheckRowSum();
    // Mdt.CheckMatrix();

    // ##########################################
    // Timestepping        (M/dt + K + C) * u_{n+1}  =  F + M/dt * u_{n}
    // ##########################################


    // ################################## SIMULATION (i) ##################################

    // Initialize temperature u_0
    vector<double> uv(SK.Nrows(), 0.0);                                                               // temperature
    mesh_c.Init_Solution_mult(uv, 0, [u0_mug](double x, double y) -> double { return u0_mug; });      // mug
    mesh_c.Init_Solution_mult(uv, 1, [u0_fluid](double x, double y) -> double { return u0_fluid; });  // fluid
    mesh_c.Init_Solution_mult(uv, 2, [u0_air](double x, double y) -> double { return u0_air; });      // air
    
    //mesh_c.Visualize(uv);

    auto t3 = system_clock::now(); // start timer
    double average_cup_temperature = u0_mug;
    double percentage_temp_reached = 0.0;
    double goal_temp = 67.0;
    double goal_perc = 60.0;

    double time_count = 0;
    while (average_cup_temperature < goal_temp)
    // while (percentage_temp_reached < goal_perc)
    //for (int step = 0; step < steps; ++step)
    {
        vector<double> G(Mdt.Nrows(), 0.0);        
        Mdt.Mult(G, uv);                                           // G = M/dt * u_{n}
        
        vector<double> H = fv;   
        for (size_t i = 0; i < Mdt.Nrows(); ++i)
        {
            H[i] += G[i];                                         // H = F + G
        }

        JacobiSolve(SK, H, uv);          // solve: (M/dt + K + C) * u_{n+1}  =  F + M/dt * u_{n}
                                         //        ----- SK -----               ------ H -------

        average_cup_temperature = mesh_c.AverageVectorFunction_perSubdomain(uv, 0);
        percentage_temp_reached = mesh_c.CheckTemp_mult(uv, 0, goal_temp);
        cout << "Average cup temperature: " << average_cup_temperature << " after " << time_count << " seconds. " << endl;
        cout << "% of elements reached temperature " << goal_temp << "ºC: " << percentage_temp_reached << endl;
        time_count += dt;
    }
    auto t4 = system_clock::now();  // stop timer


    auto duration = duration_cast<microseconds>(t4 - t3);        // duration in microseconds
    double t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds
    cout << "\n\nJacobiSolve: timing in sec. : " << t_diff << endl;

    mesh_c.Visualize(uv);

    // ################################## SIMULATION (iii) ##################################
    double u0_coffee = 85.0;
    mesh_c.Init_Solution_mult(uv, 1, [u0_coffee](double x, double y) -> double { return u0_coffee; });  // fluid

    t3 = system_clock::now(); // start timer
    double average_coffee_temperature = u0_coffee;
    percentage_temp_reached = 0.0;

    // ------------------------ initialize preCICE ------------------------
    int commRank = 0;
    int commSize = 1;
    std::string configFileName("precise-config.xml");
    std::string solverName("Solid");
    std::string meshName("Solid-Mesh");

    std::cout << "Running Solid-solver with preCICE config file \"" << configFileName << "\" and participant name \"" << solverName << "\".\n";

    Participant participantSolid(solverName, configFileName, commRank, commSize);

    int meshDim  = participantSolid.getMeshDimensions(meshName); // gets mesh dimensions (=2) from config file
    int numberOfVertices; // number of vertices at "wet" surface
    {
        // TODO: DETERMINE NUMBER OF VERTICES AT TOP
    }
    std::vector<double> coords(numberOfVertices*meshDim); // coords of vertices at "wet" surface
    {
        // TODO: DETERMINE COORDINATES OF VERTICES AT TOP
    }
    std::vector<int> vertexIDs(numberOfVertices);
    
    participantSolid.setMeshVertices(meshName, coords, vertexIDs);

    // ----- initialize read- and write-data
    std::vector<double> temperature(numberOfVertices * meshDim);    // read-data
    std::vector<double> heatFlux(numberOfVertices * meshDim);       // write-data


    double solverDt = 1.0; // solver time step size
    double preciceDt = 1.0; // maximum precice time step size
    dt = 1.0; // actual time step size

    participantSolid.initialize();


    // ------------------------ timestepping ------------------------
    goal_temp = 50.0;
    goal_perc = 60.0;

    time_count = 0;
    while (average_cup_temperature > goal_temp)
    // while (percentage_temp_reached > goal_perc)
    {
        preciceDt = participantSolid.getMaxTimeStepSize();
        solverDt = 1.0;
        dt = min(preciceDt, solverDt);

        // ----- read temperature computed by openFOAM simulation -----
        participantSolid.readData("Solid-Mesh", "Temperature", vertexIDs, dt, temperature);
        {
            // TODO: SET THE MESH TO THE INCOMING TEMPERATURE
            // "initNewTemp(temperature)"  
        }

        // ----- solve time step -----
        vector<double> G(Mdt.Nrows(), 0.0);        
        Mdt.Mult(G, uv);                                           // G = M/dt * u_{n}
        
        vector<double> H = fv;   
        for (size_t i = 0; i < Mdt.Nrows(); ++i)
        {
            H[i] += G[i];                                         // H = F + G
        }

        JacobiSolve(SK, H, uv);          // solve: (M/dt + K + C) * u_{n+1}  =  F + M/dt * u_{n}
                                         //        ----- SK -----               ------ H -------

        average_coffee_temperature = mesh_c.AverageVectorFunction_perSubdomain(uv, 1);
        percentage_temp_reached = mesh_c.CheckTemp_mult(uv, 0, goal_temp);
        cout << "Average coffee temperature: " << average_coffee_temperature << " after " << time_count << " seconds. " << endl;
        cout << "% of elements reached temperature " << goal_temp << "ºC: " << percentage_temp_reached << endl;

        // ----- write the heat-flux, so openFOAM can read it
        {
            // TODO: COMPUTE THE HEAT-FLUX FROM THE NEW uv TEMPERATURE SOLUTION
            // "computeHeatFlux(heatFlux)"
        }
        participantSolid.writeData("Solid-Mesh", "Heat-Flux", vertexIDs, heatFlux);

        // ----- advance preCICE -----
        participantSolid.advance(dt);


        time_count += dt;
    }
    t4 = system_clock::now();  // stop timer
    duration = duration_cast<microseconds>(t4 - t3);        // duration in microseconds
    t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds
    cout << "\n\nJacobiSolve: timing in sec. : " << t_diff << endl;


    // ----- free data, close communication -----
    participantSolid.finalize();



    return 0;
}