integration.cc 2.49 KB
Newer Older
Peter Bastian's avatar
Peter Bastian committed
1 2 3 4
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

// Dune includes
5
#include <config.h>             // file generated by CMake
6 7
#include <iomanip>
#include <iostream>
8
#include <dune/grid/yaspgrid.hh> // load yaspgrid definition
9
#include <dune/common/parallel/mpihelper.hh> // include mpi helper class
Peter Bastian's avatar
Peter Bastian committed
10

11 12
#include "functors.hh"
#include "integrateentity.hh"
Peter Bastian's avatar
Peter Bastian committed
13

14 15 16 17 18 19 20
//! uniform refinement test
template<class Grid>
void uniformintegration (Grid& grid)
{
  // function to integrate
  Exp<typename Grid::ctype,Grid::dimension> f;

21 22 23 24
  // get GridView on leaf grid - type
  typedef typename Grid :: LeafGridView GridView;

  // get GridView instance
25
  GridView gridView = grid.leafGridView();
26

27
  // get iterator type
28
  typedef typename GridView :: template Codim<0> :: Iterator LeafIterator;
29 30 31

  // loop over grid sequence
  double oldvalue=1E100;
32
  for (int k=0; k<10; k++)
33
  {
34 35
    // compute integral with some order
    double value = 0.0;
36 37
    LeafIterator eendit = gridView.template end<0>();
    for (LeafIterator it = gridView.template begin<0>(); it!=eendit; ++it)
38
      value += integrateEntity(*it,f,1);                /*@\label{ic:call}@*/
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54

    // print result and error estimate
    std::cout << "elements="
              << std::setw(8) << std::right
              << grid.size(0)
              << " integral="
              << std::scientific << std::setprecision(12)
              << value
              << " error=" << std::abs(value-oldvalue)
              << std::endl;

    // save value of integral
    oldvalue=value;

    // refine all elements
    grid.globalRefine(1);
Peter Bastian's avatar
Peter Bastian committed
55
  }
56
}
Peter Bastian's avatar
Peter Bastian committed
57

58
int main(int argc, char **argv)
Peter Bastian's avatar
Peter Bastian committed
59
{
Robert Klöfkorn's avatar
Robert Klöfkorn committed
60 61
  // initialize MPI, finalize is done automatically on exit
  Dune::MPIHelper::instance(argc,argv);
62

63 64
  // start try/catch block to get error messages from dune
  try {
Robert Klöfkorn's avatar
Robert Klöfkorn committed
65 66
    using namespace Dune;

67 68 69 70
    // the GridSelector :: GridType is defined in gridtype.hh and is
    // set during compilation
    typedef GridSelector :: GridType Grid;

Robert Klöfkorn's avatar
Robert Klöfkorn committed
71 72
    // use unitcube from grids
    std::stringstream dgfFileName;
73 74
    dgfFileName << DUNE_GRID_HOWTO_EXAMPLE_GRIDS_PATH
      << "unitcube" << Grid::dimension << ".dgf";
Robert Klöfkorn's avatar
Robert Klöfkorn committed
75

76 77
    // create grid pointer
    GridPtr<Grid> gridPtr( dgfFileName.str() );
78

79
    // integrate and compute error with extrapolation
Robert Klöfkorn's avatar
Robert Klöfkorn committed
80
    uniformintegration( *gridPtr );
81 82
  }
  catch (std::exception & e) {
83
    std::cout << "ERROR: " << e.what() << std::endl;
84 85 86 87 88 89
    return 1;
  }
  catch (...) {
    std::cout << "Unknown ERROR" << std::endl;
    return 1;
  }
90

Peter Bastian's avatar
Peter Bastian committed
91 92 93
  // done
  return 0;
}