Add new example

src/05-poisson-problem.cc

+// -*- 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>
+#include <dune/common/stringutility.hh>
+
+// included dune-geometry headers
+#include <dune/geometry/quadraturerules.hh>
+
+// included dune-grid headers
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#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>
+
+
+
+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;
+}
+
+
+
+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;
+
+//    using Grid = Dune::YaspGrid<2>;
+//    auto gridPointer = createCubeGrid<Grid>();
+
+//    using Grid = Dune::YaspGrid<3>;
+//    auto gridPointer = createCubeGrid<Grid>();
+
+//    using Grid = Dune::UGGrid<2>;
+//    auto gridPointer = createCubeGrid<Grid>();
+
+    using Grid = Dune::UGGrid<2>;
+    auto gridPointer = createSimplexGrid<Grid>();
+
+    Grid& grid = *gridPointer;
+
+    grid.globalRefine(5);
+    for(int i=0; i<2; ++i)
+    {
+      auto gridView = grid.levelGridView(grid.maxLevel());
+      for(const auto& e : Dune::elements(gridView))
+        if (e.geometry().center().two_norm() < .5)
+          grid.mark(1, e);
+      grid.adapt();
+    }
+    auto gridView = grid.leafGridView();
+    using GridView = decltype(gridView);
+
+
+    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double,1>>;
+
+    Matrix A;
+    Vector rhs;
+    Vector x;
+
+
+    auto rhsFunction = [](auto x) {
+      return x.two_norm() < .5;
+    };
+    assemblePoissonProblem(gridView, A, rhs, rhsFunction);
+
+    x.resize(rhs.size(), 0);
+
+    std::vector<bool> isBoundary;
+    isBoundary.resize(rhs.size(), false);
+    computeBoundaryVertices(gridView, isBoundary);
+
+    for(std::size_t i=0; i<rhs.size(); ++i)
+    {
+      if (isBoundary[i])
+      {
+        for(auto& entry: A[i])
+          entry = 0;
+        A[i][i] = 1;
+        rhs[i] = x[i];
+      }
+    }
+
+    Dune::MatrixAdapter<Matrix,Vector,Vector> op(A);
+    Dune::SeqILU0<Matrix,Vector,Vector> ilu0(A,1.0);
+    Dune::SeqSSOR<Matrix,Vector,Vector> preconditioner(A,1,1.0);
+    Dune::CGSolver<Vector> cg(op,
+        preconditioner, // preconditione
+        1e-4, // desired residual reduction factor
+        200,   // maximum number of iterations
+        2);   // verbosity of the solver
+
+    // Object storing some statistics about the solving process
+    Dune::InverseOperatorResult statistics;
+
+    // Solve!
+    cg.apply(x, rhs, statistics);
+
+    Dune::VTKWriter<GridView> vtkWriter(gridView);
+    vtkWriter.addVertexData(x, "solution");
+    vtkWriter.write("05-poisson-solution");
+
+    return 0;
+  }
+  catch (Dune::Exception &e){
+    std::cerr << "Dune reported error: " << e << std::endl;
+  }
+  catch (...){
+    std::cerr << "Unknown exception thrown!" << std::endl;
+  }
+}

src/CMakeLists.txt

@@ -12,3 +12,6 @@ target_link_dune_default_libraries("03-function-integration")
 
 add_executable("04-gridviews" 04-gridviews.cc)
 target_link_dune_default_libraries("04-gridviews")
+
+add_executable("05-poisson-problem" 05-poisson-problem.cc)
+target_link_dune_default_libraries("05-poisson-problem")