recipe-grid-function-operations.cc

See explanation at Operations on grid functions

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
 
#include <iostream>
 
#include <dune/common/filledarray.hh>
#include <dune/common/parallel/mpihelper.hh>
 
#include <dune/grid/yaspgrid.hh>
 
#include <dune/pdelab.hh>
 
int main(int argc, char** argv)
{
  // Initialize Mpi
  Dune::MPIHelper::instance(argc, argv);
 
  // need a grid in order to test grid functions
  constexpr unsigned int dim = 2;
  Dune::FieldVector<double,dim> L(5.0);
  std::array<int,dim> N(Dune::filledArray<dim,int>(64));
 
  typedef Dune::YaspGrid<dim> Grid;
  Grid grid(L,N);
 
  // [Defining functions]
  auto analyticFunction1 = Dune::PDELab::makeGridFunctionFromCallable (grid.leafGridView(), [&](const auto& x){
    return exp(-(x*x));
  });
  auto analyticFunction2 = Dune::PDELab::makeGridFunctionFromCallable (grid.leafGridView(), [&](const auto& x){
    return x[0]*3.0;
  });
 
 
  {
    // [Instantiating DifferenceSquaredAdapter]
    Dune::PDELab::DifferenceSquaredAdapter<decltype(analyticFunction1),decltype(analyticFunction2)> difference(analyticFunction1, analyticFunction2);
    // Note: Newer compiler versions do not require template arguments here!
    auto integral = Dune::PDELab::integrateGridFunction(difference,10);
    std::cout << "Integral: " << integral << std::endl;
  }
 
  {
    // [Defining operations]
    auto difference = Dune::PDELab::makeGridFunctionFromCallable (grid.leafGridView(), [&](const auto& element, const auto& xlocal){
      Dune::FieldVector<double,1> y1;
      analyticFunction1.evaluate(element, xlocal, y1);
      Dune::FieldVector<double,1> y2;
      analyticFunction2.evaluate(element, xlocal, y2);
      y2 -= y1;
      return y2.two_norm2();
    });
    // [Compute integral]
    auto integral = Dune::PDELab::integrateGridFunction(difference,10);
    std::cout << "Integral: " << integral << std::endl;
  }
}