recipe-integration.cc

See explanation at Integrate a function on a grid

// SPDX-FileCopyrightText: Copyright (C) DUNE Project contributors, see file COPYING in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// always include the config file
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
// C++ includes
#include<math.h>
#include<iostream>
// dune-common includes
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertreeparser.hh>
#include<dune/common/timer.hh>
// dune-geometry includes
#include<dune/geometry/referenceelements.hh>
#include<dune/geometry/quadraturerules.hh>
// dune-grid includes
#include <dune/grid/yaspgrid.hh>
 
 
int main(int argc, char** argv)
{
  // Maybe initialize Mpi
  [[maybe_unused]] Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
 
  // [set up grid]
  const int dim = 4;
  using Grid = Dune::YaspGrid<dim>;
  Dune::FieldVector<double,dim> len; for (auto& l:len) l=1.0;
  std::array<int,dim> cells; for (auto& c : cells) c=5;
  Grid grid(len,cells);
 
  // [small vectors and matrices]
  Dune::FieldVector<double,4> x({1,2,3,4}); // make a vector
  auto y(x); // copy constructor
  y *= 1.0/3.0; // scaling
  [[maybe_unused]] auto s = x*y; // scalar product
  [[maybe_unused]] auto norm = x.two_norm(); // Euclidean norm
  Dune::FieldMatrix<double,4,4> A({{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}); // make a matrix
  A.mv(x,y); // matvec: y = Ax
  A.usmv(0.5,x,y); // axpy: y += 0.5*Ax
 
  // [a function to integrate]
  auto u = [](const auto& x){return std::exp(x.two_norm());};
 
  // [integration with midpoint rule]
  double integral=0.0;
  auto gv = grid.leafGridView(); // extract the grid view
  for (const auto& e : elements(gv))
    integral += u(e.geometry().center())*e.geometry().volume();
  std::cout << "integral = " << integral << std::endl;
 
  // [integration with quadrature rule]
  double integral2 = 0.0;
  using QR = Dune::QuadratureRules<Grid::ctype,dim>;
  for (const auto& e : elements(gv))
    {
      auto geo = e.geometry();
      auto quadrature = QR::rule(geo.type(),5);
      for (const auto& qp : quadrature)
        integral2 += u(geo.global(qp.position()))
          *geo.integrationElement(qp.position())*qp.weight();
    }
  std::cout << "integral2 = " << integral2 << std::endl;
 
  // [integrating a flux]
  auto f = [](const auto& x){return x;};
  double divergence=0.0;
  for (const auto& i : elements(gv)) {
    for (const auto& I : intersections(gv,i))
      if (!I.neighbor())
        {
          auto geoI = I.geometry();
          divergence += f(geoI.center())*I.centerUnitOuterNormal()*geoI.volume();
        }
  }
  std::cout << "divergence = " << divergence << std::endl;
}