diff --git a/src/03-function-integration.cc b/src/03-function-integration.cc new file mode 100644 index 0000000000000000000000000000000000000000..d32971c2f12894d292ccc2b5602974301e3057cf --- /dev/null +++ b/src/03-function-integration.cc @@ -0,0 +1,117 @@ +// -*- 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 + + + +// included standard library headers +#include <iostream> +#include <array> + +// included dune-common headers +#include <dune/common/parallel/mpihelper.hh> +#include <dune/common/exceptions.hh> +#include <dune/common/fvector.hh> + +// included dune-geometry headers +#include <dune/geometry/quadraturerules.hh> + +// included dune-grid headers +#include <dune/grid/yaspgrid.hh> +#include <dune/grid/io/file/vtk/vtkwriter.hh> + + + + + +template<class GridView, class F> +double integrate(const GridView& gridView, const F& f, std::size_t order) +{ + static const int dim = GridView::dimension; + + double value = 0.0; + for(const auto& e: Dune::elements(gridView)) + { + auto geometryType = e.type(); + + const auto& geometry = e.geometry(); + + const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order); + + for (size_t i=0; i < quad.size(); ++i) + { + const auto& quadPos = quad[i].position(); + const auto& weight = quad[i].weight(); + + auto globalQuadPos = geometry.global(quadPos); + + // get integration factor + const double integrationElement = geometry.integrationElement(quadPos); + + value += f(globalQuadPos) * weight * integrationElement; + } + } + return value; +} + + +int main(int argc, char** argv) +{ + try{ + // Maybe initialize MPI + Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv); + + // Print process rank + if(Dune::MPIHelper::isFake) + std::cout<< "This is a sequential program." << std::endl; + else + std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size() + <<" processes!"<<std::endl; + + // fix grid dimension + const int dim = 2; + + // use the YaspGrid implementation + // YaspGrid = "Yet another structured parallel grid" + using Grid = Dune::YaspGrid<dim>; + + // extension of the domain [0,l] + Dune::FieldVector<double,dim> l(1); + + // start with 10x10x... elements + std::array<int,dim> E; + for(auto& e : E) + e = 2; + + // create grid + Grid grid(l, E); + + // refine the grid globally two times + grid.globalRefine(7); + + // create function to integrate + auto f = [](const auto& x) { + return x.infinity_norm(); + }; + + // get a leaf GridView + auto gridView = grid.leafGridView(); + + std::cout << integrate(gridView, f, 7) << std::endl; + + using GridView = decltype(gridView); + Dune::VTKWriter<GridView> vtkWriter(gridView); + vtkWriter.write("01-refined-grid"); + + return 0; + } + catch (Dune::Exception &e){ + std::cerr << "Dune reported error: " << e << std::endl; + } + catch (...){ + std::cerr << "Unknown exception thrown!" << std::endl; + } +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 72a1a439c3ec615245369ea6e9c9a33c578778ba..bbeefd9a7a2412cd0480b55f126c98ecda51ca23 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,3 +6,6 @@ target_link_dune_default_libraries("01-basic-application") add_executable("02-basic-grid-interface" 02-basic-grid-interface.cc) target_link_dune_default_libraries("02-basic-grid-interface") + +add_executable("03-function-integration" 03-function-integration.cc) +target_link_dune_default_libraries("03-function-integration")