New example on global interpolation

src/06-interpolation.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>
+#include <math.h>
+
+// 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-functions headers
+#include <dune/functions/common/functionfromcallable.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 Vector, class F>
+void interpolateFunction(const GridView& gridView, Vector& v, const F& f)
+{
+  static const int dim = GridView::dimension;
+  const auto& indexSet = gridView.indexSet();
+  std::size_t size = indexSet.size(dim);
+
+  v.resize(size);
+  v = 0;
+
+  Dune::PQkLocalFiniteElementCache<double, double, dim, 1> feCache;
+
+  for(const auto& e: Dune::elements(gridView))
+  {
+    const auto& geometry = e.geometry();
+    const auto& localFiniteElement = feCache.get(e.type());
+    auto localSize = localFiniteElement.localBasis().size();
+
+    using Element = std::decay_t<decltype(e)>;
+    using LocalDomain = typename Element::Geometry::LocalCoordinate;
+    using Range = double;
+
+    auto localF = [&](const LocalDomain& x) {
+      return f(geometry.global(x));
+    };
+
+    std::vector<Range> localV;
+    
+//    localFiniteElement.localInterpolation().interpolate(localF, localV);
+
+    using FunctionFromCallable = typename Dune::Functions::FunctionFromCallable<Range(LocalDomain), decltype(localF)>;
+    auto localFWithEvaluate = FunctionFromCallable(localF);
+    localFiniteElement.localInterpolation().interpolate(localFWithEvaluate, localV);
+
+    for(std::size_t i=0; i<localSize; ++i)
+    {
+      std::size_t globalI = indexSet.subIndex(e, localFiniteElement.localCoefficients().localKey(i).subEntity(), dim);
+      v[globalI] = localV[i];
+    }
+  }
+}
+
+
+
+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)*100.0;
+    };
+    assemblePoissonProblem(gridView, A, rhs, rhsFunction);
+
+    x.resize(rhs.size(), 0);
+
+    auto dirichletFunction = [](auto x) {
+      return sin(x[0] * 2.0*M_PI)*.1;
+    };
+    interpolateFunction(gridView, x, dirichletFunction);
+
+    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("06-interpolation");
+
+    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

@@ -15,3 +15,6 @@ target_link_dune_default_libraries("04-gridviews")
 
 add_executable("05-poisson-problem" 05-poisson-problem.cc)
 target_link_dune_default_libraries("05-poisson-problem")
+
+add_executable("06-interpolation" 06-interpolation.cc)
+target_link_dune_default_libraries("06-interpolation")