#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/mapping_q.h>

#include <fstream>
#include <iostream>

using namespace dealii;

class GmeshTest
{

public:
    static constexpr int dim  = 2 ;

public:

    GmeshTest(bool _use_hypercube)
    :
        triangulation(),
        dof_handler(triangulation),
        use_hypercube(_use_hypercube)
    {
        make_grid(); 
        test_mapping();
    }

private:
    void make_grid()
    {

        std::string outmesh_name;
        if(use_hypercube)
        {
            std::vector< unsigned  int > repetitions(dim);
            repetitions[0] = 60;
            repetitions[1] = 20;

            std::array<double,dim> L = {6.0, 2.0};
            GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point<dim>(0.0,0.0), Point<dim>(L[0],L[1]), true);
            outmesh_name="mesh_hypercube.vtk";
        }
        else
        {
            GridIn<dim> gridin;
            gridin.attach_triangulation(triangulation);
            std::ifstream f("test.msh"); 
            gridin.read_msh(f);
            outmesh_name="mesh_gmsh.vtk";
        }

        std::ofstream out(outmesh_name);
        GridOut grid_out;
        grid_out.write_vtk(triangulation, out);
    
    }

    void test_mapping()
    {
        auto cell = dof_handler.begin_active();
        auto endc = dof_handler.end();
        for (; cell!=endc; ++cell)
        {
            auto cellN = cell->index();
            if (cellN == 354)
            {
                Point<dim> test;
                test(0) = 5.45;
                test(1) = 0.55;

                MappingQ<dim> mapping(1);

                Point<dim> dp = mapping.transform_real_to_unit_cell(cell,test);
                Point<dim> dp_real = mapping. transform_unit_to_real_cell(cell,dp);

                std::cout << "cell=" << cell->index() << std::endl;
                std::cout << "cell vertex(0): " << cell->vertex(0) << std::endl;
                std::cout << "cell vertex(1): " << cell->vertex(1) << std::endl;
                std::cout << "cell vertex(2): " << cell->vertex(2) << std::endl;
                std::cout << "cell vertex(3): " << cell->vertex(3) << std::endl;
                std::cout << "=================================================" << std::endl;
                std::cout << "test_point=   " << test << "  mapped_point=   " << dp <<" backmapping: "<<dp_real<< std::endl;
                std::cout << "=================================================" << std::endl;
            }
        }
    }

private:
    Triangulation<dim>      triangulation;
    DoFHandler<dim>         dof_handler;
    const bool              use_hypercube;
};

int main(int argc, char *argv[])
{

    GmeshTest hypercube(true);
    GmeshTest gmsh(false);

    return 0;
}