#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#include "moab/Core.hpp"
#include "Coupler.hpp"
#include "moab_mpi.h"
#include "ElemUtil.hpp"


using namespace moab;
using namespace std;

void createMoabMesh(Interface * mb, string tagName, Tag & vtag, EntityHandle & src_mesh, EntityHandle & tar_mesh)
{
    double src_ps[] = {-1, -1, 0,  0, -1, 0,  1, -1, 0,
                        -1, 0, 0, 0, 0, 0, 1, 0, 0,
                        -1, 1, 0, 0, 1, 0, 1, 1, 0};

    double tar_ps[] = {-1, -1, 0,  0, -1, 0,  1, -1, 0,
                        -1, 0, 0, 0.3, 0.4, 0, 1, 0, 0,
                        -1, 1, 0, 0, 1, 0, 1, 1, 0};

    //double src_data[] = {0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0};
    double src_data[] = {1.0, 2.0, 2.0, 3.0};

    Range src_vertices;
    Range tar_vertices;
    mb->create_vertices(src_ps,  9,  src_vertices);
    mb->create_vertices(tar_ps,  9,  tar_vertices);


    int cells[] = {0, 1, 4, 3, 1, 2, 5, 4, 3, 4, 7, 6, 4, 5, 8, 7};
    int cellLocations[] = {0, 4, 8, 12};

    vector<EntityHandle>  src_cells(4, 0);
    vector<EntityHandle> tar_cells(4, 0);
    for(int i =0; i < 4; i++)
    {
        int idx0 = cells[cellLocations[i]];
        int idx1 = cells[cellLocations[i] + 1];
        int idx2 = cells[cellLocations[i] + 2];
        int idx3 = cells[cellLocations[i] + 3];
        EntityHandle src_conn[] = {src_vertices[idx0], src_vertices[idx1], src_vertices[idx2], src_vertices[idx3]};
        mb->create_element(MBQUAD, src_conn, 4, src_cells[i]);

        EntityHandle tar_conn[] = {tar_vertices[idx0], tar_vertices[idx1], tar_vertices[idx2], tar_vertices[idx3]};
        mb->create_element(MBQUAD, tar_conn, 4, tar_cells[i]);
    }

    // Create tag
    double value = 0.0;
    int value_len = 1;
    mb->tag_get_handle(tagName.c_str(), value_len, MB_TYPE_DOUBLE, vtag,
                       MB_TAG_DENSE|MB_TAG_CREAT, &value);

    mb->tag_set_data(vtag, &src_cells[0], 4, src_data);

    mb->create_meshset(MESHSET_SET, src_mesh);
    mb->add_entities(src_mesh, src_vertices);
    mb->add_entities(src_mesh, src_cells.data(), src_cells.size());

    mb->create_meshset(MESHSET_SET, tar_mesh);
    mb->add_entities(tar_mesh, tar_vertices);

    mb->add_entities(tar_mesh, tar_cells.data(), tar_cells.size());
}

int main(int argc, char *argv[])
{
    int err = MPI_Init(&argc, &argv);
    Interface * mb = new Core();
    ParallelComm * pc0 = new ParallelComm(mb, MPI_COMM_WORLD);


    // prepare test mesh and test data
    EntityHandle src_mesh=0;
    EntityHandle tar_mesh=0;
    Tag vtag;
    string tagName = "rho";

    createMoabMesh(mb, tagName, vtag, src_mesh, tar_mesh);

    Range src_elems, tar_elems;
    mb->get_entities_by_dimension(src_mesh, 2, src_elems);
    mb->get_entities_by_dimension(tar_mesh, 2, tar_elems);

    cout << " Source elem size:" << src_elems.size() << endl;
    cout << " Target elem size:" << tar_elems.size() << endl;
    cout << "Test 1" << endl;
     Coupler mbc(mb, pc0, src_elems, 0, true, 2);

     ErrorCode result;
     double toler = 5.e-10;
     int nprocs, rank;
     int ierr = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
     ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);

     cout << nprocs << ", " << rank << endl;

     Coupler::Method method = Coupler::CONSTANT;

     // Get points  coordinates in targets mesh
     std::vector<double> vpos;
     int num_of_points = (int)tar_elems.size();
     vpos.resize(3*tar_elems.size());
     mb->get_coords(tar_elems, &vpos[0]);

     for(int i = 0 ; i < num_of_points; i++)
         cout << vpos[3*i]<< ", " << vpos[3*i+1] << ", " << vpos[3*i+2] << endl;

     // locate the position
     mbc.locate_points(vpos.data(), num_of_points, 0, toler);

     // Compute the interpolation value on source mesh
     std::vector<double> field(num_of_points);
     mbc.interpolate(method, vtag, &field[0]);

     // Set field values as tag on target vertices
     mb->tag_set_data(vtag, tar_elems, &field[0]);

     for(int i = 0; i < 4; i++)
     {
         cout << field[i] << endl;
     }

     delete mb;
     delete pc0;

     // Finalize the MPI environment.
     err = MPI_Finalize();

     return 0;

}