interpolation example added

dune.module

@@ -8,4 +8,4 @@ Version: 0.1
 Maintainer: simon.praetorius@tu-dresden.de
 #depending on
 Depends: dune-common dune-grid
-Suggests: dune-uggrid dune-alugrid dune-foamgrid dune-functions dune-geometry
+Suggests: dune-uggrid dune-alugrid dune-foamgrid dune-functions dune-geometry dune-istl

dune/multimesh/mmiterator.hh

@@ -176,7 +176,7 @@ namespace Dune
     MultiMeshMasterLeafIterator (tag::begin_iterator, const GridImp* multiMesh, std::size_t master)
       : Super{tag::begin_iterator{}, multiMesh}
       , incrementAllowed_(multiMesh->size(), true)
-      , level_(multiMesh->size())
+      , level_(multiMesh->size(), 0)
       , multiEntity_(multiMesh->size())
       , master_(master)
     {
@@ -225,7 +225,7 @@ namespace Dune
     // got down to leaf entity on master grid and until the master grid level for the other grids
     bool levelReached_impl (std::size_t i, const HostEntity& entity) const
     {
-      return (i == master_ && entity.isLeaf()) || entity.level() == level_[master_];
+      return entity.isLeaf() || (i != master_ && entity.level() == level_[master_]);
     }
 
   private:

src/CMakeLists.txt

@@ -16,6 +16,12 @@ target_link_dune_default_libraries("multigridview")
 add_executable("hierarchiciterator" hierarchiciterator.cc)
 target_link_dune_default_libraries("hierarchiciterator")
 
+add_executable("interpolate" interpolate.cc)
+target_link_dune_default_libraries("interpolate")
+if (HAVE_ALBERTA)
+  add_dune_alberta_flags("interpolate")
+endif ()
+
 find_package(MTL PATHS /opt/development/mtl4)
 if (MTL_FOUND)
   set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")

src/interpolate.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
+
+#if ! HAVE_DUNE_FUNCTIONS
+#error "Require dune-functions!"
+#endif
+
+#if ! HAVE_DUNE_ALUGRID
+#error "Require dune-alugrid!"
+#endif
+
+#if ! HAVE_ALBERTA
+#error "Require Alberta!"
+#endif
+
+#include <iostream>
+#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
+#include <dune/common/exceptions.hh> // We use exceptions
+#include <dune/common/timer.hh>
+
+#include <dune/alugrid/grid.hh>
+
+#include <dune/multimesh/mmhierarchiciterator.hh>
+#include <dune/multimesh/mmgridfactory.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/grid/albertagrid.hh>
+#include <dune/grid/albertagrid/albertareader.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+#include <dune/istl/bvector.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+using namespace Dune;
+using namespace Dune::Indices;
+
+using VectorType = BlockVector<FieldVector<double,1>>;
+
+template <class VectorType, class Basis>
+void writeFile(VectorType const& x, Basis const& basis, std::string filename)
+{
+  auto gv = basis.gridView();
+  auto xFct = Functions::makeDiscreteGlobalBasisFunction<double>(basis, TypeTree::hybridTreePath(), x);
+  VTKWriter<decltype(gv)> vtkWriter(gv);
+  vtkWriter.addVertexData(xFct, VTK::FieldInfo("u", VTK::FieldInfo::Type::scalar, 1));
+  vtkWriter.write(filename);
+}
+
+template <class Grids, class SignedDistFct>
+void markElements(Grids& grids, SignedDistFct&& signedDistFct, double eps = 0.1)
+{
+  for (auto const& elems : elements(grids.leafGridView())) {
+    auto geo = elems[0].geometry();
+    int level = elems[0].level();
+    int newLevel = 0;
+    for (int j = 0; j < geo.corners(); ++j)
+      newLevel = std::max(newLevel, std::abs(signedDistFct(geo.corner(j))) < eps ? 12 : 0);
+
+    int m = level < newLevel ? 1 : -1;
+    grids[0].mark(m, elems[0]);
+    grids[1].mark(m, elems[1]);
+  }
+}
+
+template <class Element, class GlobalCoordinate>
+auto findLeafEntity(Element const& father, GlobalCoordinate const& global)
+{
+  const int childLevel = father.level()+1;
+
+  // loop over all child Entities
+  const auto end = father.hend(childLevel);
+  for (auto it = father.hbegin(childLevel); it != end; ++it)
+  {
+    Element child = *it;
+    auto geo = child.geometry();
+
+    auto local = geo.local(global);
+    if (referenceElement(geo).checkInside(local))
+    {
+      // return if we found the leaf, else search through the child entites
+      if (child.isLeaf())
+        return std::make_pair(std::move(child), local);
+      else
+        return findLeafEntity(child, global);
+    }
+  }
+  assert(father.isLeaf());
+  return std::make_pair(father, father.geometry().local(global));
+}
+
+template <class Grids, class OldBasis, class NewBasis>
+void mmInterpolate(Grids const& grids,
+                   VectorType const& oldCoeff, OldBasis const& oldBasis,
+                   VectorType& newCoeff, NewBasis const& newBasis)
+{
+  using namespace Dune::Indices;
+  using LocalCoordinate = typename Grids::Traits::template Codim<0>::Geometry::LocalCoordinate;
+  auto oldLocalView = oldBasis.localView();
+  auto newLocalView = newBasis.localView();
+
+  newCoeff.resize(newBasis.dimension());
+  for (const auto& e : master_leaf_elements(grids, 1))
+  {
+    auto const& oldEntity = e[0];
+    auto const& newEntity = e[1]; // always a leaf entity
+    assert(newEntity.isLeaf());
+
+    newLocalView.bind(newEntity);
+    auto newGeo = newEntity.geometry();
+    auto const& node = newLocalView.tree();
+    auto const& fe = node.finiteElement();
+
+    std::size_t newFeSize = fe.size();
+    std::vector<double> newDOFs(newFeSize);
+
+    // entity not changed
+    if (newEntity.level() == oldEntity.level() && oldEntity.isLeaf()) {
+      oldLocalView.bind(oldEntity);
+      std::size_t oldFeSize = oldLocalView.tree().finiteElement().size();
+      assert(newFeSize == oldFeSize);
+      for (std::size_t i = 0; i < newFeSize; ++i)
+        newDOFs[i] = oldCoeff[oldLocalView.index(oldLocalView.tree().localIndex(i))];
+    }
+    // entity was coarsened
+    if (newEntity.level() == oldEntity.level() && !oldEntity.isLeaf()) {
+      auto evalLeaf = [&](LocalCoordinate const& x) {
+        thread_local std::vector<FieldVector<double,1>> shapeValues;
+        auto leafLocal = findLeafEntity(oldEntity, newGeo.global(x));
+
+        oldLocalView.bind(leafLocal.first);
+        auto const& node = oldLocalView.tree();
+        auto const& fe = node.finiteElement();
+        fe.localBasis().evaluateFunction(leafLocal.second, shapeValues);
+
+        double y = 0;
+        for (std::size_t i = 0; i < shapeValues.size(); ++i)
+          y += shapeValues[i] * oldCoeff[oldLocalView.index(node.localIndex(i))];
+
+        return y;
+      };
+
+      auto const& node = newLocalView.tree();
+      auto const& fe = node.finiteElement();
+
+      using FFC = Functions::FunctionFromCallable<double(LocalCoordinate), decltype(evalLeaf)>;
+      fe.localInterpolation().interpolate(FFC(evalLeaf), newDOFs);
+    }
+    // entity was refined
+    else if (newEntity.level() >= oldEntity.level()) {
+      oldLocalView.bind(oldEntity);
+      auto oldGeo = oldEntity.geometry();
+      auto evalOld = [&](LocalCoordinate const& x)
+      {
+        thread_local std::vector<FieldVector<double,1>> shapeValues;
+        auto const& fe = oldLocalView.tree().finiteElement();
+        fe.localBasis().evaluateFunction(oldGeo.local(newGeo.global(x)), shapeValues);
+
+        double y = 0;
+        for (std::size_t i = 0; i < shapeValues.size(); ++i)
+          y += shapeValues[i] * oldCoeff[oldLocalView.index(node.localIndex(i))];
+
+        return y;
+      };
+
+      using FFC = Functions::FunctionFromCallable<double(LocalCoordinate), decltype(evalOld)>;
+      fe.localInterpolation().interpolate(FFC(evalOld), newDOFs);
+    }
+
+    for (std::size_t i = 0; i < newFeSize; ++i)
+      newCoeff[newLocalView.index(node.localIndex(i))] = newDOFs[i];
+  }
+}
+
+int main(int argc, char** argv)
+{
+  MPIHelper::instance(argc, argv);
+
+  assert(argc > 1);
+  std::string filename = argv[1];
+
+  // using HostGrid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
+  using HostGrid = Dune::AlbertaGrid<2,2>;
+  using WorldVector = typename HostGrid::template Codim<0>::Geometry::GlobalCoordinate;
+  using Grid = MultiMesh<HostGrid>;
+  AlbertaReader<Grid> albertaReader;
+  GridFactory<Grid> factory(2);
+  albertaReader.readGrid(filename, factory);
+
+  std::unique_ptr<Grid> gridPtr(factory.createGrid());
+  auto& grids = *gridPtr;
+  grids.globalRefine(6);
+
+  double eps = 0.1;
+  WorldVector shift(0.0);
+  auto signedDistFct = [&shift](WorldVector const& x) { return (x - shift).two_norm() - 0.7; };
+
+  // initial refinement
+  for (int i = 0; i < 8; ++i) {
+    markElements(grids, signedDistFct, eps);
+    grids.preAdapt();
+    grids.adapt();
+    grids.postAdapt();
+  }
+
+  using namespace Functions::BasisBuilder;
+  auto basis0 = makeBasis(grids.leafGridView(0), lagrange<2>()); // old grid
+  auto basis1 = makeBasis(grids.leafGridView(1), lagrange<2>()); // new grid
+
+  // interpolate phase-field to fine grid
+  VectorType oldPhase(basis0.dimension());
+  VectorType phase(basis1.dimension());
+  Functions::interpolate(basis1, phase, [eps,d=signedDistFct](auto const& x)
+  {
+    return 0.5*(1.0 - std::tanh(2.0*d(x)/eps));
+  });
+  writeFile(phase, basis1, "phase_0");
+
+  Timer t;
+  double timeMarkElements = 0;
+  double timeAdapt = 0;
+  double timeBasisUpdate = 0;
+  double timeInterpolate = 0;
+  double timeWriteFile = 0;
+
+  // assemble a sequence of systems for finer and finer grids
+  for (int i = 1; i < 20; ++i) {
+    std::cout << i << "\n";
+    oldPhase.resize(basis0.dimension());
+    oldPhase = phase;
+
+    shift[0] += 0.01;
+
+    t.reset();
+    markElements(grids, signedDistFct, eps);
+    timeMarkElements += t.elapsed();
+
+    t.reset();
+    grids[1].preAdapt();
+    grids[1].adapt();
+    grids[1].postAdapt();
+    timeAdapt += t.elapsed();
+
+    t.reset();
+    basis1.update(grids.leafGridView(1));
+    timeBasisUpdate += t.elapsed();
+
+    t.reset();
+    mmInterpolate(grids, oldPhase, basis0, phase, basis1); // interpolate from old grid to new grid
+    timeInterpolate += t.elapsed();
+
+    t.reset();
+    writeFile(phase, basis1, "phase_" + std::to_string(i));
+    timeWriteFile += t.elapsed();
+
+    t.reset();
+    grids[0].preAdapt();
+    grids[0].adapt();
+    grids[0].postAdapt();
+    timeAdapt += t.elapsed();
+
+    t.reset();
+    basis0.update(grids.leafGridView(0));
+    timeBasisUpdate += t.elapsed();
+  }
+
+  std::cout << "timings:\n";
+  std::cout << "  time(markElements) = " << timeMarkElements << "\n";
+  std::cout << "  time(adapt) = " << timeAdapt << "\n";
+  std::cout << "  time(basisUpdate) = " << timeBasisUpdate << "\n";
+  std::cout << "  time(interpolate) = " << timeInterpolate << "\n";
+  std::cout << "  time(writeFile) = " << timeWriteFile << "\n";
+}