Add new example

src/07-poisson-problem-functions.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/io/file/vtk/subsamplingvtkwriter.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-functions headers
+#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+// included dune-fu-tutorial headers
+#include <dune/fu-tutorial/referenceelementutility.hh>
+
+
+// included dune-fu-tutorial headers
+#include "04-gridviews.hh"
+#include "05-poisson-problem.hh"
+
+
+template<class GridView, class RHSFunction, class Basis>
+void assemblePoissonProblem(
+    const GridView& gridView,
+    Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>& matrix,
+    Dune::BlockVector<Dune::FieldVector<double,1>>& rhs,
+    const RHSFunction& rhsFunction,
+    const Basis& basis)
+{
+  std::size_t size = basis.size();
+
+  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, 100, 0.1);
+
+  ElementMatrix elementMatrix;
+  ElementRhs elementRhs;
+
+  auto localView = basis.localView();
+  auto localIndexSet = basis.localIndexSet();
+
+  for(const auto& e: Dune::elements(gridView))
+  {
+    localView.bind(e);
+    localIndexSet.bind(localView);
+
+    const auto& localFiniteElement = localView.tree().finiteElement();
+
+    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 = localIndexSet.index(i);
+      for(std::size_t j=0; j<localSize; ++j)
+      {
+        std::size_t globalJ = localIndexSet.index(j);
+        matrixBuilder[globalI][globalJ] += elementMatrix[i][j];
+      }
+
+      rhs[globalI] += elementRhs[i];
+    }
+  }
+  matrix.compress();
+}
+
+template<class GridView, class BitVector, class Basis>
+void computeBoundaryVertices(const GridView& gridView, BitVector& isBoundary, const Basis& basis)
+{
+  using namespace Dune;
+  using namespace Dune::FuTutorial;
+
+  auto localView = basis.localView();
+  auto localIndexSet = basis.localIndexSet();
+
+  for(const auto& element: elements(gridView))
+  {
+    localView.bind(element);
+    localIndexSet.bind(localView);
+
+    const auto& localFiniteElement = localView.tree().finiteElement();
+    std::size_t localSize = localFiniteElement.localBasis().size();
+
+    for(const auto& intersection: intersections(gridView, element))
+    {
+      if (intersection.boundary())
+      {
+        for(std::size_t i=0; i<localSize; ++i)
+        {
+          auto localKey = localFiniteElement.localCoefficients().localKey(i);
+          if (localKey.codim() > 0)
+            for(const auto& subEntity: subEntities(referenceElement(element), insideFacet(intersection), codim(localKey.codim())))
+              if (subEntity == localKey.subEntity())
+                isBoundary[localIndexSet.index(i)] = 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(1);
+    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)*10.0;
+    };
+
+
+    using Basis = typename Dune::Functions::PQkNodalBasis<GridView,4>;
+    Basis basis(gridView);
+
+    assemblePoissonProblem(gridView, A, rhs, rhsFunction, basis);
+
+    x.resize(basis.size(), 0);
+
+    auto dirichletFunction = [](auto x) {
+      return sin(x[0] * 2.0*M_PI)*.1;
+    };
+//    Dune::Functions::interpolate(basis, Dune::TypeTree::hybridTreePath(), dirichletNodes, dirichletFunction);
+    Dune::Functions::interpolate(basis, x, dirichletFunction);
+
+    std::vector<bool> isBoundary;
+    isBoundary.resize(basis.size(), false);
+    computeBoundaryVertices(gridView, isBoundary, basis);
+
+    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);
+
+    auto xFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(basis, x);
+
+   {
+      Dune::VTKWriter<GridView> vtkWriter(gridView);
+      vtkWriter.addVertexData(xFunction, Dune::VTK::FieldInfo("solution", Dune::VTK::FieldInfo::Type::scalar, 1));
+      vtkWriter.write("07-poisson-solution-functions");
+    } {
+      Dune::SubsamplingVTKWriter<GridView> vtkWriter(gridView,4);
+      vtkWriter.addVertexData(xFunction, Dune::VTK::FieldInfo("solution", Dune::VTK::FieldInfo::Type::scalar, 1));
+      vtkWriter.write("07-poisson-solution-functions-subsamples");
+    }
+    
+
+    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

@@ -18,3 +18,6 @@ target_link_dune_default_libraries("05-poisson-problem")
 
 add_executable("06-interpolation" 06-interpolation.cc)
 target_link_dune_default_libraries("06-interpolation")
+
+add_executable("07-poisson-problem-functions" 07-poisson-problem-functions.cc)
+target_link_dune_default_libraries("07-poisson-problem-functions")