#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include <iostream>
#include <vector>
#include <array>

#include <dune/common/exceptions.hh>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/grid/uggrid.hh>             // Koristimo UGGrid
#include <dune/grid/io/file/gmshreader.hh> // GmshReader klasa
#include <dune/grid/io/file/vtk.hh>
#include <dune/geometry/quadraturerules.hh>

template <int dim>
Dune::FieldVector<double, dim> q(Dune::FieldVector<double, dim>const & x)
{
    auto y = x;
    y /= x.dot(x);
    return y;
}

template<typename GV>
double integration(GV gv, int p)
{
    const int dim = 2;
    double integral = 0.0;
    for (auto const & element : elements(gv))
    {
        for (auto const & side : intersections(gv, element))
        {
            if (side.boundary())
            {
                const auto side_geom = side.geometry();
                auto outer_normal = side.centerUnitOuterNormal();
                const auto &rule = Dune::QuadratureRules<double, dim - 1>::rule(side_geom.type(), p);

                for (auto const &qpoint : rule)
                {
                    auto fval = q(side_geom.global(qpoint.position()));;
                    double weight = qpoint.weight();
                    double detjac = side_geom.integrationElement(qpoint.position());
                    integral += (fval * outer_normal) * weight * detjac;
                }
            }
        }
    }

    return integral;
}

int main(int argc, char** argv)
{
    const int dim = 2;

    Dune::MPIHelper::instance(argc, argv);
    using GridType = Dune::UGGrid<dim>;
    using GridView = GridType::LeafGridView;

    if (argc < 2)
        std::cout << "Need msh file name!" << std::endl, std::exit(-1);

    GridType *grid = Dune::GmshReader<GridType>::read(argv[1]);
    grid->loadBalance();

    GridView gv = grid->leafGridView();

    auto &set = gv.indexSet();
    std::vector<double> q_data(dim * gv.size(dim));

    for (auto const & vertex : vertices(gv))
        for (unsigned int i = 0; i < dim; ++i)
            q_data[dim * set.index(vertex) + i] = q(vertex.geometry().center())[i];

    Dune::VTKWriter<GridView> writer(gv);
    writer.addVertexData(q_data, "q", dim);
    writer.write("flux");

    std::cout << "Flux = " << integration(gv, 5) << std::endl;
    return 0;
}