//		MPI code in C++.
//		See [Gropp/Lusk/Skjellum, "Using MPI", p.33/41 etc.]
//		and  /opt/mpich/include/mpi2c++/comm.h  for details

#include "elements.h"
#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>
using namespace std;
using namespace std::chrono;  // timing

int main(int /* argc */, char ** /* argv */ )
{
//#define PDE_2D
#ifdef PDE_2D
    //Mesh const mesh("square_100.txt");
    Mesh const mesh("square_tiny.txt");
    //mesh.Debug();
    //mesh.DebugEdgeBased();
    //constexpr P1_2d_1dof elem;
    //const P1_2d_1dof elem;

    const P1_2d elem;
    FEM_Matrix_2 SK(mesh, elem);           // CRS matrix

    vector<double> uv(SK.Nrows(), 0.0);    // temperature
    vector<double> fv(SK.Nrows(), 0.0);    // r.h.s.        
    //SK.CalculateLaplace(fv);
    SK.CalculateLaplace(fv, rhs_lap2);
    //SK.CheckMatrix();
    //SK.Debug();
    //mesh.SetValues(uv, [](double x, double y) -> double {
    //return x * x * std::sin(2.5 * M_PI * y);
    //} );
    mesh.SetValues(uv, sol_lap2);

    SK.ApplyDirichletBC(uv, fv); // Dirichlet nodes are defines in input file

    //////SK.writeBinary("sparseMatrix.bin");
    //SK.Debug();
#else
    //Mesh mesh("NS_3D.txt");
    Mesh mesh("GH_NS_3D.txt");
    FEM_Matrix_2 SK(mesh, P1_3d());                  // CRS matrix
    //mesh.liftToQuadratic();
    ////FEM_Matrix_2 SK(mesh,P2_3d(3));                 // CRS matrix
    //FEM_Matrix_2 SK(mesh,P1_2vec_3d());               // CRS matrix

    vector<double> uv(SK.Nrows(), 0.0);    // temperature
    vector<double> fv(SK.Nrows(), 0.0);    // r.h.s.
    //SK.CalculateLaplace(fv);
    SK.CalculateLaplace(fv,rhs_lap3);
    //SK.CheckMatrix();
    //mesh.SetValues(uv, [](double x, double y, double z) -> double {
    //return x * x * std::sin(2.5 * M_PI * y) + z*z;
    //} );
    mesh.SetValues(uv, sol_lap3);

    SK.ApplyDirichletBC_Box(uv, fv, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0); // brute force Dirichlet via box
    //SK.ApplyDirichletBC_Box(uv,fv,5.0,6.0,5.0,6.0,0.0,0.0); // z==0 ==> Dirichlet-BC
#endif


    //SK.Debug();
    auto exact_sol(uv);

    ////SK.Mult(fv,uv);

    //auto t3 = system_clock::now(); // start timer
    JacobiSolve(SK, fv, uv );          // solve the system of equations
    //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 << "JacobiSolve: timing in sec. : " << t_diff << endl;

    auto [val, idx] = findLargestAbsError(exact_sol, uv, 1e-6, 10);
    cout << val << " :idx: " << idx << endl;

    //cout << "exact: " << exact_sol << endl;
    //cout << "numer: " << uv << endl;

    ////mesh.Visualize(getAbsError(exact_sol, uv));


    ////mesh.Write_ascii_matlab("uv.txt", uv);
    mesh.Visualize(uv);

    return 0;
}