// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#include <functional>
#include <iostream>
#include <numeric>

#include <dune/common/filledarray.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/test/testsuite.hh>

#include <dune/grid/yaspgrid.hh>

#include <dune/multimesh/multimesh.hh>

using namespace Dune;

template <std::size_t dim, class Test>
bool test_dim(Test& test)
{
  FieldVector<double,dim> lower; lower = -1.5;
  FieldVector<double,dim> upper; upper =  1.5;
  auto num_elements = filledArray<dim>(2);

  using HostGrid = YaspGrid<dim, EquidistantOffsetCoordinates<double,dim>>;
  MultiMesh<HostGrid> grid(3, lower, upper, num_elements);

  // volume of the domain
  double domain = std::inner_product(lower.begin(), lower.end(), upper.begin(), 1.0, std::multiplies<>{},
    [](double l, double u) { return std::abs(u - l); });
  std::cout << "volume(domain<" << dim << ">) = " << domain << "\n";

  grid[0].globalRefine(2);
  grid[1].globalRefine(1);
  grid[2].globalRefine(3);

  // calculate volume by summing up the entity volumes of the smalles leaf entities
  double volume_leaf = 0.0;
  for (auto const& entities : elements(grid.leafGridView())) {
    auto it_small = std::max_element(entities.begin(), entities.end(),
      [](auto const& e1, auto const& e2) { return e1.level() < e2.level(); });

    auto geo_small = it_small->geometry();
    volume_leaf += geo_small.volume();
  }
  std::cout << "volume_leaf(elements<" << dim << ">) = " << volume_leaf << "\n";
  test.check(std::abs(volume_leaf - domain) < 1.e-10);

  // calculate volume by summing up the entity volumes of the level=1 entities
  double volume_level = 0.0;
  for (auto const& entities : elements(grid.levelGridView(1))) {
    auto geo = entities[0].geometry();
    volume_level += geo.volume();
  }
  std::cout << "volume_level(elements<" << dim << ">) = " << volume_level << "\n";
  test.check(std::abs(volume_level - domain) < 1.e-10);
}

int main(int argc, char** argv)
{
  MPIHelper::instance(argc, argv);

  Dune::TestSuite test;
  test_dim<1>(test);
  test_dim<2>(test);
  test_dim<3>(test);

  return test.exit();
}