#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 RefineTest
{

public:
    static constexpr int dim  = 2 ;

public:

    RefineTest(bool _refine)
    :
        triangulation(),
        dof_handler(triangulation),
        refine(_refine)
    {
        make_grid(); 
        test_mapping();
    }

private:
    void make_grid()
    {
        if(refine)
        {
            std::vector< unsigned  int > repetitions(dim);
            repetitions[0] = repetitions[1] = 5;
            std::array<double,dim> L = {1.0, 1.0};
            GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point<dim>(0.0,0.0), Point<dim>(L[0],L[1]), true);
 	    triangulation.refine_global(1);
        }
        else
        {
            std::vector< unsigned  int > repetitions(dim);
            repetitions[0] = repetitions[1] = 10;
            std::array<double,dim> L = {1.0, 1.0};
            GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point<dim>(0.0,0.0), Point<dim>(L[0],L[1]), true);
 	    triangulation.refine_global(0); 
        }
    }

    void test_mapping()
    {
        auto cell = dof_handler.begin_active();
        auto endc = dof_handler.end();
        for (; cell!=endc; ++cell)
        {

            auto cellN = cell->index();
            auto P1 = cell->vertex(0);
	    auto P2 = cell->vertex(1);
	    auto P3 = cell->vertex(2);

	    Point<dim> test;
            test(0) = (P1[0]+P2[0])/2.0; 
            test(1) = (P2[1]+P3[1])/2.0; 

            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);

	    if (abs(test[0]-dp_real[0]) > 1e-8 || abs(test[1]-dp_real[1]) > 1e-8)
	    {
		std::cout << " " << std::endl;
		std::cout << "ERROR" << std::endl;
                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              refine;
};

int main(int argc, char *argv[])
{
    RefineTest refineglobal(true);
    return 0;
}