From 70153cc753b16cce4ed2102b4bfcd3ee569eb025 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@mi.fu-berlin.de>
Date: Thu, 1 Jun 2017 09:35:35 +0200
Subject: [PATCH] New example
---
src/03-function-integration.cc | 117 +++++++++++++++++++++++++++++++++
src/CMakeLists.txt | 3 +
2 files changed, 120 insertions(+)
create mode 100644 src/03-function-integration.cc
diff --git a/src/03-function-integration.cc b/src/03-function-integration.cc
new file mode 100644
index 0000000..d32971c
--- /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 72a1a43..bbeefd9 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")
--
GitLab