diff --git a/src/05-poisson-problem.cc b/src/05-poisson-problem.cc
index 23baf461bc04a1657844c90375df02b0ff362e22..0dc6d8cb39e95df528284511b2db978426ed257b 100644
--- a/src/05-poisson-problem.cc
+++ b/src/05-poisson-problem.cc
@@ -41,200 +41,9 @@
#include <dune/fu-tutorial/referenceelementutility.hh>
-
-template<class Grid>
-auto createCubeGrid()
-{
- static const int dim = Grid::dimension;
-
- auto lowerLeft = Dune::FieldVector<double,dim>(0);
- auto upperRight = Dune::FieldVector<double,dim>(1);
- auto elementsPerDirection = std::array<unsigned int,dim>();
- for(auto& e : elementsPerDirection)
- e = 2;
-
- return Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, elementsPerDirection);
-}
-
-template<class Grid>
-auto createSimplexGrid()
-{
- static const int dim = Grid::dimension;
-
- auto lowerLeft = Dune::FieldVector<double,dim>(0);
- auto upperRight = Dune::FieldVector<double,dim>(1);
- auto elementsPerDirection = std::array<unsigned int,dim>();
- for(auto& e : elementsPerDirection)
- e = 2;
-
- return Dune::StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, elementsPerDirection);
-}
-
-
-template<class GridView>
-void writeGridView(const GridView& gridView, std::string postFix)
-{
- Dune::VTKWriter<GridView> vtkWriter(gridView);
- vtkWriter.write(std::string("04-gridviews-")+postFix);
-}
-
-
-
-// Compute local per element matrix
-template <class Element, class LocalFiniteElement, class MatrixType>
-void assembleLocalStiffnessMatrix(const Element& element, const LocalFiniteElement& localFiniteElement, MatrixType& elementMatrix)
-{
- auto geometryType = element.type();
- const auto& geometry = element.geometry();
- const int dim = Element::dimension;
- std::size_t localSize = localFiniteElement.localBasis().size();
-
- // Resize and initialize matrix
- elementMatrix.setSize(localSize, localSize);
- elementMatrix = 0; // fills the entire matrix with zeroes
-
- // Get a quadrature rule
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order);
-
- // Loop over all quadrature points
- for (size_t k=0; k < quad.size(); ++k)
- {
- const auto& quadPos = quad[k].position();
- const auto& weight = quad[k].weight();
-
- // The transposed inverse Jacobian of the map from the reference element to the element
- const auto& jacobian = geometry.jacobianInverseTransposed(quadPos);
-
- // The multiplicative factor in the integral transformation formula
- const auto& integrationElement = geometry.integrationElement(quadPos);
-
- // The gradients of the shape functions on the reference element
- std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
-
- // Compute the shape function gradients on the real basis functions
- std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());
- for (size_t i=0; i<gradients.size(); i++)
- jacobian.mv(referenceGradients[i][0], gradients[i]);
-
- // Compute the actual matrix entries
- for (size_t i=0; i<elementMatrix.N(); i++)
- for (size_t j=0; j<elementMatrix.M(); j++ )
- elementMatrix[i][j] += ( gradients[i] * gradients[j] ) * weight * integrationElement;
- }
-}
-
-
-
-// Compute local per element right hand side
-template <class Element, class LocalFiniteElement, class LocalRHSVector, class RHSFunction>
-void assembleLocalRHS(
- const Element& element,
- const LocalFiniteElement& localFiniteElement,
- LocalRHSVector& localRhs,
- const RHSFunction& rhsFunction)
-{
- auto geometryType = element.type();
- const auto& geometry = element.geometry();
- const int dim = Element::dimension;
- std::size_t localSize = localFiniteElement.localBasis().size();
-
- // Resize and initialize vector
- localRhs.resize(localSize);
- localRhs = 0;
-
- // Get a quadrature rule
- int order = 2*(dim*localFiniteElement.localBasis().order());
- const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order);
-
- // Loop over all quadrature points
- for (size_t k=0; k < quad.size(); ++k)
- {
- const auto& quadPos = quad[k].position();
- const auto& weight = quad[k].weight();
-
- // Quadrature position in global world coordinates
- const auto& globalQuadPos = geometry.global(quadPos);
-
- // The multiplicative factor in the integral transformation formula
- const auto& integrationElement = geometry.integrationElement(quadPos);
-
- double functionValue = rhsFunction(globalQuadPos);
-
- // Evaluate all shape function values at this point
- std::vector<Dune::FieldVector<double,1> > shapeFunctionValues;
- localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
-
- // Actually compute the vector entries
- for (size_t i=0; i<localSize; i++)
- localRhs[i] += shapeFunctionValues[i] * functionValue * weight * integrationElement;
- }
-}
-
-
-template<class GridView, class RHSFunction>
-void assemblePoissonProblem(
- const GridView& gridView,
- Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>& matrix,
- Dune::BlockVector<Dune::FieldVector<double,1>>& rhs,
- const RHSFunction& rhsFunction)
-{
- static const int dim = GridView::dimension;
- const auto& indexSet = gridView.indexSet();
- std::size_t size = indexSet.size(dim);
-
- using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
- using ElementMatrix = Dune::Matrix<Dune::FieldMatrix<double,1,1>>;
- using ElementRhs = Dune::BlockVector<Dune::FieldVector<double,1>>;
-
- rhs.resize(size);
- rhs = 0;
- Dune::ImplicitMatrixBuilder<Matrix> matrixBuilder(matrix, size, size, 10, 0.1);
-
- Dune::PQkLocalFiniteElementCache<double, double, dim, 1> feCache;
-
- ElementMatrix elementMatrix;
- ElementRhs elementRhs;
-
- for(const auto& e: Dune::elements(gridView))
- {
- const auto& localFiniteElement = feCache.get(e.type());
- std::size_t localSize = localFiniteElement.localBasis().size();
-
- assembleLocalStiffnessMatrix(e, localFiniteElement, elementMatrix);
- assembleLocalRHS(e, localFiniteElement, elementRhs, rhsFunction);
-
- for(std::size_t i=0; i<localSize; ++i)
- {
- std::size_t globalI = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(i).subEntity(), dim);
- for(std::size_t j=0; j<localSize; ++j)
- {
- std::size_t globalJ = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(j).subEntity(), dim);
- matrixBuilder[globalI][globalJ] += elementMatrix[i][j];
- }
-
- rhs[globalI] = elementRhs[i];
- }
- }
- matrix.compress();
-}
-
-
-
-template<class GridView, class BitVector>
-void computeBoundaryVertices(const GridView& gridView, BitVector& isBoundary)
-{
- using namespace Dune;
- using namespace Dune::FuTutorial;
-
- const auto& indexSet = gridView.indexSet();
- for(const auto& element: elements(gridView))
- for(const auto& intersection: intersections(gridView, element))
- if (intersection.boundary())
- for(const auto& vertex: subVertices(referenceElement(element), insideFacet(intersection)))
- isBoundary[subIndex(indexSet, element, vertex)] = true;
-}
+// included dune-fu-tutorial headers
+#include "04-gridviews.hh"
+#include "05-poisson-problem.hh"
diff --git a/src/05-poisson-problem.hh b/src/05-poisson-problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5d8d6129157d65ab1a4de7c76a43eda9ac1b16b2
--- /dev/null
+++ b/src/05-poisson-problem.hh
@@ -0,0 +1,197 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FUTUTORIAL_SRC_05POISSON_PROBLEM_HH
+#define DUNE_FUTUTORIAL_SRC_05POISSON_PROBLEM_HH
+
+// 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>
+#include <dune/common/stringutility.hh>
+
+// included dune-geometry headers
+#include <dune/geometry/quadraturerules.hh>
+
+// included dune-grid headers
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/uggrid.hh>
+
+// included dune-istl headers
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+
+
+// included dune-localfunctions headers
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+// included dune-fu-tutorial headers
+#include <dune/fu-tutorial/referenceelementutility.hh>
+
+
+
+// Compute local per element matrix
+template <class Element, class LocalFiniteElement, class MatrixType>
+void assembleLocalStiffnessMatrix(const Element& element, const LocalFiniteElement& localFiniteElement, MatrixType& elementMatrix)
+{
+ auto geometryType = element.type();
+ const auto& geometry = element.geometry();
+ const int dim = Element::dimension;
+ std::size_t localSize = localFiniteElement.localBasis().size();
+
+ // Resize and initialize matrix
+ elementMatrix.setSize(localSize, localSize);
+ elementMatrix = 0; // fills the entire matrix with zeroes
+
+ // Get a quadrature rule
+ int order = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order);
+
+ // Loop over all quadrature points
+ for (size_t k=0; k < quad.size(); ++k)
+ {
+ const auto& quadPos = quad[k].position();
+ const auto& weight = quad[k].weight();
+
+ // The transposed inverse Jacobian of the map from the reference element to the element
+ const auto& jacobian = geometry.jacobianInverseTransposed(quadPos);
+
+ // The multiplicative factor in the integral transformation formula
+ const auto& integrationElement = geometry.integrationElement(quadPos);
+
+ // The gradients of the shape functions on the reference element
+ std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+ // Compute the shape function gradients on the real basis functions
+ std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());
+ for (size_t i=0; i<gradients.size(); i++)
+ jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+ // Compute the actual matrix entries
+ for (size_t i=0; i<elementMatrix.N(); i++)
+ for (size_t j=0; j<elementMatrix.M(); j++ )
+ elementMatrix[i][j] += ( gradients[i] * gradients[j] ) * weight * integrationElement;
+ }
+}
+
+
+
+// Compute local per element right hand side
+template <class Element, class LocalFiniteElement, class LocalRHSVector, class RHSFunction>
+void assembleLocalRHS(
+ const Element& element,
+ const LocalFiniteElement& localFiniteElement,
+ LocalRHSVector& localRhs,
+ const RHSFunction& rhsFunction)
+{
+ auto geometryType = element.type();
+ const auto& geometry = element.geometry();
+ const int dim = Element::dimension;
+ std::size_t localSize = localFiniteElement.localBasis().size();
+
+ // Resize and initialize vector
+ localRhs.resize(localSize);
+ localRhs = 0;
+
+ // Get a quadrature rule
+ int order = 2*(dim*localFiniteElement.localBasis().order());
+ const auto& quad = Dune::template QuadratureRules<double, dim>::rule(geometryType, order);
+
+ // Loop over all quadrature points
+ for (size_t k=0; k < quad.size(); ++k)
+ {
+ const auto& quadPos = quad[k].position();
+ const auto& weight = quad[k].weight();
+
+ // Quadrature position in global world coordinates
+ const auto& globalQuadPos = geometry.global(quadPos);
+
+ // The multiplicative factor in the integral transformation formula
+ const auto& integrationElement = geometry.integrationElement(quadPos);
+
+ double functionValue = rhsFunction(globalQuadPos);
+
+ // Evaluate all shape function values at this point
+ std::vector<Dune::FieldVector<double,1> > shapeFunctionValues;
+ localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
+
+ // Actually compute the vector entries
+ for (size_t i=0; i<localSize; i++)
+ localRhs[i] += shapeFunctionValues[i] * functionValue * weight * integrationElement;
+ }
+}
+
+
+template<class GridView, class RHSFunction>
+void assemblePoissonProblem(
+ const GridView& gridView,
+ Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>& matrix,
+ Dune::BlockVector<Dune::FieldVector<double,1>>& rhs,
+ const RHSFunction& rhsFunction)
+{
+ static const int dim = GridView::dimension;
+ const auto& indexSet = gridView.indexSet();
+ std::size_t size = indexSet.size(dim);
+
+ using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
+ using ElementMatrix = Dune::Matrix<Dune::FieldMatrix<double,1,1>>;
+ using ElementRhs = Dune::BlockVector<Dune::FieldVector<double,1>>;
+
+ rhs.resize(size);
+ rhs = 0;
+ Dune::ImplicitMatrixBuilder<Matrix> matrixBuilder(matrix, size, size, 10, 0.1);
+
+ Dune::PQkLocalFiniteElementCache<double, double, dim, 1> feCache;
+
+ ElementMatrix elementMatrix;
+ ElementRhs elementRhs;
+
+ for(const auto& e: Dune::elements(gridView))
+ {
+ const auto& localFiniteElement = feCache.get(e.type());
+ std::size_t localSize = localFiniteElement.localBasis().size();
+
+ assembleLocalStiffnessMatrix(e, localFiniteElement, elementMatrix);
+ assembleLocalRHS(e, localFiniteElement, elementRhs, rhsFunction);
+
+ for(std::size_t i=0; i<localSize; ++i)
+ {
+ std::size_t globalI = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(i).subEntity(), dim);
+ for(std::size_t j=0; j<localSize; ++j)
+ {
+ std::size_t globalJ = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(j).subEntity(), dim);
+ matrixBuilder[globalI][globalJ] += elementMatrix[i][j];
+ }
+
+ rhs[globalI] = elementRhs[i];
+ }
+ }
+ matrix.compress();
+}
+
+
+
+template<class GridView, class BitVector>
+void computeBoundaryVertices(const GridView& gridView, BitVector& isBoundary)
+{
+ using namespace Dune;
+ using namespace Dune::FuTutorial;
+
+ const auto& indexSet = gridView.indexSet();
+ for(const auto& element: elements(gridView))
+ for(const auto& intersection: intersections(gridView, element))
+ if (intersection.boundary())
+ for(const auto& vertex: subVertices(referenceElement(element), insideFacet(intersection)))
+ isBoundary[subIndex(indexSet, element, vertex)] = true;
+}
+
+
+
+#endif // DUNE_FUTUTORIAL_SRC_05POISSON_PROBLEM_HH