testvolume.cc 2.18 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
// -*- 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>
14
#include <dune/common/test/testsuite.hh>
15 16 17 18 19 20 21

#include <dune/grid/yaspgrid.hh>

#include <dune/multimesh/multimesh.hh>

using namespace Dune;

22 23
template <std::size_t dim, class Test>
bool test_dim(Test& test)
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
{
  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
42
  double volume_leaf = 0.0;
43 44 45 46 47
  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();
48
    volume_leaf += geo_small.volume();
49
  }
50 51
  std::cout << "volume_leaf(elements<" << dim << ">) = " << volume_leaf << "\n";
  test.check(std::abs(volume_leaf - domain) < 1.e-10);
52

53 54 55 56 57 58 59 60
  // 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);
61 62 63 64 65
}

int main(int argc, char** argv)
{
  MPIHelper::instance(argc, argv);
66 67 68 69 70 71 72

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

  return test.exit();
73
}