MultiMeshHierarchicIterator added to traverse leaf-childs of an element

dune/multimesh/mmhierarchiciterator.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_MULTIMESH_HIERARCHIC_ITERATOR_HH
+#define DUNE_MULTIMESH_HIERARCHIC_ITERATOR_HH
+
+#include <numeric>
+#include <stack>
+#include <variant>
+
+#include <dune/common/iteratorfacades.hh>
+#include <dune/common/std/type_traits.hh>
+
+namespace Dune
+{
+  /** \brief Iterator over all entities of a given codimension and level of a grid.
+   *  \ingroup MultiMesh
+   */
+  template <class HostGrid>
+  class MultiMeshHierarchicIterator
+      : public ForwardIteratorFacade<MultiMeshHierarchicIterator<HostGrid>,
+            typename HostGrid::template Codim<0>::Entity,
+            typename HostGrid::template Codim<0>::Entity>
+  {
+  private:
+    template <int codim, PartitionIteratorType pitype, class HostGrid_>
+    friend class MultiMeshIterator;
+
+    using HostEntity = typename HostGrid::template Codim<0>::Entity;
+    using EntityTest = std::function<bool(HostEntity)>;
+
+    struct EntityStackEntry
+    {
+      template <class Entity>
+      explicit EntityStackEntry (Entity&& entity)
+        : it(entity.hbegin(entity.level()+1))
+        , end(entity.hend(entity.level()+1))
+      {}
+
+      typename HostEntity::HierarchicIterator it;
+      typename HostEntity::HierarchicIterator end;
+
+      friend bool operator== (const EntityStackEntry& lhs, const EntityStackEntry& rhs)
+      {
+        return lhs.it == rhs.it;
+      }
+    };
+
+    class EntityStack
+        : public std::stack<EntityStackEntry, std::vector<EntityStackEntry>>
+    {
+      using Super = std::stack<EntityStackEntry, std::vector<EntityStackEntry>>;
+
+    public:
+      explicit EntityStack (int maxLevel)
+      {
+        c.reserve(maxLevel);
+      }
+
+      // return true if all levels >= l are finished, i.e. hierarchic iterators it+1 == end
+      bool finished (std::size_t l = 0) const
+      {
+        bool f = true;
+        for (auto tmp = c[l].it; f && l < c.size(); ++l) {
+          tmp = c[l].it;
+          f = f && (++tmp) == c[l].end;
+        }
+        return f;
+      }
+
+      friend bool operator== (const EntityStack& lhs, const EntityStack& rhs)
+      {
+        return lhs.size() == rhs.size() && lhs.c == rhs.c;
+      }
+
+    protected:
+      using Super::c;
+    };
+
+  public:
+    /// Constructor. Stores a pointer to the entity
+    template <class Entity, class GridView>
+    MultiMeshHierarchicIterator (const Entity& entity, const GridView& gridView)
+      : entity_{&entity}
+      , contains_{[gv=gridView](const HostEntity& e) { return gv.contains(e); }}
+      , maxLevel_{gridView.grid().maxLevel()}
+      , entityStack_{EntityStack{gridView.grid().maxLevel()}}
+    {
+      initialIncrement();
+    }
+
+    /// Constructor which create the end iterator
+    /**
+     * \param  endDummy   Here only to distinguish it from the other constructor
+     */
+    MultiMeshHierarchicIterator (bool endDummy)
+      : entityStack_{0}
+    {}
+
+    /// prefix increment
+    void increment ()
+    {
+      // 1. go up in tree or to next entity on current level until we can go down again
+      while (!entityStack_.empty()) {
+        auto& top = entityStack_.top();
+        ++top.it;
+        if (top.it == top.end) {
+          entityStack_.pop();
+        } else {
+          break;
+        }
+      }
+
+      // 2. if entityStack is empty, iteration is finished
+      if (entityStack_.empty())
+        return;
+
+      // 3. go down in tree until leaf entity
+      auto child = dereference();
+      for (; !(contains_(child) || child.isLeaf()); child = dereference()) {
+        assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
+        entityStack_.emplace(child);
+        assert( entityStack_.size() <= maxLevel_ );
+      }
+
+      assert(contains_(child) && "No valid child element found in gridView");
+    }
+
+    /// dereferencing
+    decltype(auto) dereference () const
+    {
+      if (entityStack_.empty()) {
+        return *entity_;
+      } else {
+        assert(entityStack_.top().it != entityStack_.top().end);
+        return *entityStack_.top().it;
+      }
+    }
+
+    /// equality
+    bool equals (const MultiMeshHierarchicIterator& that) const
+    {
+      return entityStack_ == that.entityStack_;
+    }
+
+  protected:
+
+    // got to first leaf entity on gridView
+    void initialIncrement ()
+    {
+      // 1. go down in tree until leaf entity
+      auto child = dereference();
+      for (; !(contains_(child) || child.isLeaf()); child = dereference()) {
+        assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
+        entityStack_.emplace(child);
+        assert( entityStack_.size() <= maxLevel_ );
+      }
+
+      assert(contains_(child) && "No valid child element found in gridView");
+    }
+
+  private:
+    const HostEntity* entity_ = nullptr;
+
+    EntityTest contains_;
+    int maxLevel_;
+    EntityStack entityStack_;
+  };
+
+  template <class Entity, class GridView>
+  inline auto childs(const Entity& entity, const GridView& gridView)
+  {
+    using Iterator = MultiMeshHierarchicIterator<typename GridView::Grid>;
+    return IteratorRange<Iterator>{ Iterator{entity, gridView}, Iterator{true} };
+  }
+
+}  // end namespace Dune
+
+#endif // DUNE_MULTIMESH_HIERARCHIC_ITERATOR_HH

dune/multimesh/mmiterator.hh

@@ -13,14 +13,10 @@
 #include <dune/grid/common/gridenums.hh>
 
 #include "mmentity.hh"
+#include "mmhierarchiciterator.hh"
 
 namespace Dune
 {
-  namespace tag
-  {
-    struct with_gridview {};
-  }
-
   /** \brief Iterator over all entities of a given codimension and level of a grid.
    *  \ingroup MultiMesh
    */
@@ -55,45 +51,8 @@ namespace Dune
 
     using HostEntity = typename HostGrid::template Codim<0>::Entity;
 
-    using EntityTest = std::function<bool(HostEntity)>;
-
-    struct EntityStackEntry
-    {
-      template <class Entity>
-      explicit EntityStackEntry (Entity&& entity)
-        : it(entity.hbegin(entity.level()+1))
-        , end(entity.hend(entity.level()+1))
-      {}
-
-      typename HostEntity::HierarchicIterator it;
-      typename HostEntity::HierarchicIterator end;
-    };
-
-    class EntityStack
-        : public std::stack<EntityStackEntry, std::vector<EntityStackEntry>>
-    {
-      using Super = std::stack<EntityStackEntry, std::vector<EntityStackEntry>>;
-
-    public:
-      explicit EntityStack (int maxLevel)
-      {
-        c.reserve(maxLevel);
-      }
-
-      // return true if all levels >= l are finished, i.e. hierarchic iterators it+1 == end
-      bool finished (std::size_t l = 0) const
-      {
-        bool f = true;
-        for (auto tmp = c[l].it; f && l < c.size(); ++l) {
-          tmp = c[l].it;
-          f = f && (++tmp) == c[l].end;
-        }
-        return f;
-      }
-
-    protected:
-      using Super::c;
-    };
+    using EntityTest  = typename MultiMeshHierarchicIterator<HostGrid>::EntityTest;
+    using EntityStack = typename MultiMeshHierarchicIterator<HostGrid>::EntityStack;
 
   public:
     /// Constructor. Stores a pointer to the grid
@@ -142,9 +101,10 @@ namespace Dune
       : MultiMeshIterator{multiMesh, -1, endDummy}
     {}
 
+    /// Construct an iterator from n gridViews of possibly different type
     template <class... GridViews,
       std::enable_if_t<not Std::disjunction<std::is_same<GridViews,bool>...>::value, int> = 0>
-    MultiMeshIterator (tag::with_gridview, GridViews const&... gridViews)
+    MultiMeshIterator (GridViews const&... gridViews)
       : incrementAllowed_(sizeof...(GridViews), true)
       , maxLevel_{gridViews.grid().maxLevel()...}
       , multiEntity_(sizeof...(GridViews))
@@ -162,11 +122,13 @@ namespace Dune
     }
 
     template <class... GridViews>
-    MultiMeshIterator (tag::with_gridview, bool endDummy, GridViews const&... gridViews)
+    MultiMeshIterator (bool endDummy, GridViews const&... gridViews)
       : macroIterators_{gridViews.grid().levelGridView(0).template end<0,pitype>()...}
     {}
 
 #if DUNE_HAVE_CXX_VARIANT
+    /// Construct an iterator from a vector of GridViews, where each GridView can be
+    /// either a leaf- or a level-gridView, thus we need to use a variant.
     template <class... GV>
     MultiMeshIterator (const std::vector<std::variant<GV...>>& gridViews)
       : incrementAllowed_(gridViews.size(), true)
@@ -258,11 +220,14 @@ namespace Dune
       }
 
       // 3. go down in tree until leaf entity
-      for (auto child = dereference(i); !entityTest(i,child);
-                child = dereference(i)) {
+      auto child = dereference(i);
+      for (; !(contains_[i](child) || child.isLeaf()); child = dereference(i)) {
+        assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
         entityStack.emplace(child);
         assert( entityStack.size() <= maxLevel_[i] );
       }
+
+      assert(contains_[i](child) && "No valid child element found in gridView");
     }
 
     /// Return true, if all stacks with size > stack[i].size are finished
@@ -282,16 +247,19 @@ namespace Dune
       auto& macroIt = macroIterators_[i];
       auto const& macroEnd = macroEndIterators_[i];
 
-      assert( entityStack.empty() );
+      assert(entityStack.empty());
       if (macroIt == macroEnd)
         return;
 
-      // 2. go down in tree until leaf entity
-      for (auto child = dereference(i); !entityTest(i,child);
-                child = dereference(i)) {
+      // 1. go down in tree until leaf entity
+      auto child = dereference(i);
+      for (; !(contains_[i](child) || child.isLeaf()); child = dereference(i)) {
+        assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
         entityStack.emplace(child);
-        assert( entityStack.size() <= maxLevel_[i] );
+        assert(entityStack.size() <= maxLevel_[i]);
       }
+
+      assert(contains_[i](child) && "No valid child element found in gridView");
     }
 
     HostEntity dereference (std::size_t i) const
@@ -305,11 +273,6 @@ namespace Dune
       }
     }
 
-    bool entityTest (std::size_t i, HostEntity const& entity) const
-    {
-      return contains_[i](entity) || entity.isLeaf() || !entity.isRegular();
-    }
-
   private:
     std::vector<bool> incrementAllowed_;
     std::vector<EntityTest> contains_;
@@ -322,16 +285,14 @@ namespace Dune
     std::vector<EntityStack> entityStacks_;
   };
 
-  template <class> class MultiMesh;
-
   template <class... GridViews>
   inline auto multi_elements(GridViews const&... gridViews)
   {
     using GridView0 = std::tuple_element_t<0,std::tuple<GridViews...>>;
     using Iterator = MultiMeshIterator<0,All_Partition,typename GridView0::Grid>;
     using Range = IteratorRange<Iterator>;
-    return Range{ Iterator{tag::with_gridview{}, gridViews...},
-                  Iterator{tag::with_gridview{}, true, gridViews...} };
+    return Range{ Iterator{gridViews...},
+                  Iterator{true, gridViews...} };
   }
 
 }  // end namespace Dune

dune/multimesh/utility/multibasis.hh

@@ -16,7 +16,7 @@ namespace Dune
   template <class MultiMesh, class... Bases>
   struct MultiBasis
   {
-    using MultiGridView = typename MultiMesh::LeafGridView;
+    using MultiLeafGridView = typename MultiMesh::LeafGridView;
 
     /// Constructor. Stores a leafgridView of the MultiMesh and LocalView/LocalIndexSet of the bases
     MultiBasis (MultiMesh const& mm, Bases&&... bases)
@@ -27,7 +27,7 @@ namespace Dune
     {}
 
     /// Return the leafGridView of the MultiMesh
-    const MultiGridView& gridView () const { return multiGridView_; }
+    const MultiLeafGridView& gridView () const { return multiGridView_; }
 
     /// Return the I'th basis
     template <std::size_t I>
@@ -63,7 +63,7 @@ namespace Dune
     }
 
   private:
-    MultiGridView multiGridView_;
+    MultiLeafGridView multiGridView_;
     std::tuple<Bases...> bases_;
 
     std::tuple<typename Bases::LocalView...> localViews_;

src/CMakeLists.txt

@@ -13,6 +13,9 @@ target_link_dune_default_libraries("uggrid")
 add_executable("multigridview" multigridview.cc)
 target_link_dune_default_libraries("multigridview")
 
+add_executable("hierarchiciterator" hierarchiciterator.cc)
+target_link_dune_default_libraries("hierarchiciterator")
+
 find_package(MTL PATHS /opt/development/mtl4)
 if (MTL_FOUND)
   set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
@@ -25,12 +28,18 @@ if (MTL_FOUND)
     list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
   endif ()
 
-  add_executable("phasefield" phasefield.cc)
-  target_link_dune_default_libraries("phasefield")
-  if (HAVE_ALBERTA)
-    add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "phasefield")
-  endif (HAVE_ALBERTA)
-  target_include_directories("phasefield" PRIVATE ${MTL_INCLUDE_DIRS})
-  target_compile_definitions("phasefield" PRIVATE ${MTL_COMPILE_DEFINITIONS})
-  target_compile_options("phasefield" PRIVATE -Wno-deprecated-declarations)
+  set(MTL_TARGETS "")
+  list(APPEND MTL_TARGETS "phasefield" "phasefield2" "phasefield3" "phasefield4")
+
+  foreach(target ${MTL_TARGETS})
+    add_executable(${target} ${target}.cc)
+    target_link_dune_default_libraries(${target})
+    if (HAVE_ALBERTA)
+      add_dune_alberta_flags(${target})
+    endif (HAVE_ALBERTA)
+    target_include_directories(${target} PRIVATE ${MTL_INCLUDE_DIRS})
+    target_compile_definitions(${target} PRIVATE ${MTL_COMPILE_DEFINITIONS})
+    target_compile_options(${target} PRIVATE -Wno-deprecated-declarations)
+  endforeach()
+
 endif ()
\ No newline at end of file

src/hierarchiciterator.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
+
+#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>
+
+#if HAVE_DUNE_ALUGRID
+#include <dune/alugrid/grid.hh>
+#endif
+
+#include <dune/multimesh/mmhierarchiciterator.hh>
+#include <dune/multimesh/multigridview.hh>
+#include <dune/multimesh/utility/structuredgridbuilder.hh>
+
+using namespace Dune;
+
+template <class CoarseGridView, class FineGridView>
+void printGrid(const CoarseGridView& coarse, const FineGridView& fine)
+{
+  volatile std::size_t n = 0;
+  Dune::Timer t;
+  for (auto const& coarse_entity : elements(coarse)) {
+    for (auto const& fine_entity : childs(coarse_entity, fine)) {
+      n += fine.indexSet().index(fine_entity);
+      std::cout << "indices = [" << coarse.indexSet().index(coarse_entity) << ", "
+                                 << fine.indexSet().index(fine_entity) << "]\n";
+    }
+  }
+  std::cout << n << "\n";
+  std::cout << "time: " << t.elapsed() << "\n";
+}
+
+
+int main(int argc, char** argv)
+{
+  MPIHelper::instance(argc, argv);
+
+#if HAVE_DUNE_ALUGRID
+  FieldVector<double,2> lower_left = {-1.5, -1.5};
+  FieldVector<double,2> bbox = {1.5, 1.5};
+  std::array<unsigned int,2> num_elements = {2, 2};
+  using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
+
+  using Factory = StructuredGridBuilder<Grid>;
+  auto gridPtr = Factory::createSimplexGrid(lower_left, bbox, num_elements);
+  auto& grid = *gridPtr;
+
+  grid.globalRefine(4);
+
+  printGrid(grid.levelGridView(2), grid.leafGridView());
+
+#endif
+}

src/phasefield2.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/alugrid/grid.hh>
+
+#include <dune/multimesh/mmhierarchiciterator.hh>
+#include <dune/multimesh/mmgridfactory.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/grid/albertagrid/albertareader.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+#include <boost/numeric/mtl/mtl.hpp>
+#include <boost/numeric/mtl/interface/umfpack_solve.hpp>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include "phasefield2.hh"
+
+using namespace Dune;
+using namespace Dune::Indices;
+
+struct MTLVector
+{
+  using value_type = double;
+
+  MTLVector(mtl::dense_vector<double>& vec)
+    : vec_(vec) {}
+
+  void resize(std::size_t s) { vec_.change_dim(s); }
+
+  std::size_t size() const { return mtl::size(vec_); }
+
+  decltype(auto) operator[](std::size_t i) { return vec_[i]; }
+  decltype(auto) operator[](std::size_t i) const { return vec_[i]; }
+
+  decltype(auto) operator[](ReservedVector<long unsigned int, 1> i) { return vec_[i[0]]; }
+  decltype(auto) operator[](ReservedVector<long unsigned int, 1> i) const { return vec_[i[0]]; }
+
+  mtl::dense_vector<double>& vec_;
+};
+
+template <class VectorType, class Basis>
+void writeFile(VectorType& x, Basis const& basis, std::string filename)
+{
+  auto gv = basis.gridView();
+  auto xWrapper = MTLVector(x);
+  auto xFct = Functions::makeDiscreteGlobalBasisFunction<double>(basis, TypeTree::hybridTreePath(), xWrapper);
+  VTKWriter<decltype(gv)> vtkWriter(gv);
+  vtkWriter.addVertexData(xFct, VTK::FieldInfo("u", VTK::FieldInfo::Type::scalar, 1));
+  vtkWriter.write(filename);
+}
+
+
+int main(int argc, char** argv)
+{
+  MPIHelper::instance(argc, argv);
+
+  assert(argc > 1);
+  std::string filename = argv[1];
+
+  using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
+  AlbertaReader<Grid> albertaReader;
+  GridFactory<Grid> factory;
+  albertaReader.readGrid(filename, factory);
+
+  std::unique_ptr<Grid> gridPtr(factory.createGrid());
+  auto& grid = *gridPtr;
+
+  double eps = 0.05;
+  auto signedDistFct = [](auto const& x) { return x.two_norm() - 1.0; };
+
+  grid.globalRefine(10);
+  auto maxLevel = grid.maxLevel();
+
+  // refine grid[1] along phase-field interface
+  for (int i = 0; i < 4; ++i) {
+    std::cout << "refine " << i << "...\n";
+    for (auto const& elem : elements(grid.leafGridView())) {
+      auto geo = elem.geometry();
+      bool refine = false;
+      for (int j = 0; j < geo.corners(); ++j)
+        refine = refine || std::abs(signedDistFct(geo.corner(j))) < std::max(1,4-i)*eps;
+
+      if (refine)
+        grid.mark(1, elem);
+    }
+
+    grid.preAdapt();
+    grid.adapt();
+    grid.postAdapt();
+  }
+
+  using namespace Functions::BasisBuilder;
+  auto fine_basis = makeBasis(grid.leafGridView(), lagrange<1>());
+
+  using VectorType = mtl::dense_vector<double>;
+  using MatrixType = mtl::compressed2D<double>;
+
+  // interpolate phase-field to fine grid
+  VectorType phase(fine_basis.dimension());
+  auto phaseWrapper = MTLVector(phase);
+  interpolate(fine_basis, phaseWrapper, [eps,d=signedDistFct](auto const& x)
+  {
+    return 0.5*(1.0 - std::tanh(3.0*d(x)/eps));
+  });
+  writeFile(phase, fine_basis, "phase");
+
+  // assemble a sequence of systems for finer and finer grids
+  for (int l = 2; l < maxLevel; ++l) {
+    auto coarse_basis = makeBasis(grid.levelGridView(l), lagrange<1>());
+
+    VectorType rhs(coarse_basis.dimension());
+    MatrixType matrix(coarse_basis.dimension(), coarse_basis.dimension());
+    assembleSystem(coarse_basis, fine_basis, phase, matrix, rhs, eps);
+
+    VectorType u(coarse_basis.dimension());
+    u = 0.0;
+    umfpack_solve(matrix, u, rhs);
+
+    writeFile(u, coarse_basis, "u_" + std::to_string(l));
+  }
+}

src/phasefield2.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#pragma once
+
+#include <vector>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <boost/numeric/mtl/mtl.hpp>
+
+using namespace Dune;
+
+template <class LocalView, class FineBasis>
+void getLocalSystem(const LocalView& localView, const FineBasis& fine_basis,
+                    mtl::dense_vector<double> const& phase,
+                    mtl::dense2D<double>& elementMatrix,
+                    mtl::dense_vector<double>& elementVector,
+                    double eps)
+{
+  using namespace Dune::Indices;
+
+  using MultiElement = MultiEntity<typename LocalView::GridView::Grid>;
+  const int dim = MultiElement::dimension;
+  elementMatrix.change_dim(localView.size(), localView.size());
+  elementVector.change_dim(localView.size());
+  set_to_zero(elementMatrix); // fills the entire matrix with zeroes
+  elementVector = 0;
+
+  // geometry of the basis element
+  auto geometry = localView.element().geometry();
+  const auto& localFE = localView.tree().finiteElement();
+
+  auto phaseLocalView = fine_basis.localView();
+  auto phaseIndexSet = fine_basis.localIndexSet();
+  for (auto const& e : childs(localView.element(), fine_basis.gridView()))
+  {
+    phaseLocalView.bind(e);
+    phaseIndexSet.bind(phaseLocalView);
+
+    const auto& phaseLocalFE = phaseLocalView.tree().finiteElement();
+
+    auto quadGeometry = e.geometry();
+    const auto& quad = QuadratureRules<double, dim>::rule(quadGeometry.type(), 2*localFE.localBasis().order());
+    for (const auto& quadPoint : quad)
+    {
+      auto local = MultiElement::local(e, localView.element(), quadPoint.position());
+
+      const auto jacobian = geometry.jacobianInverseTransposed(local);
+      const double dx = quadGeometry.integrationElement(quadPoint.position());
+
+      static std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+      localFE.localBasis().evaluateJacobian(local, referenceGradients);
+
+      static std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+      for (std::size_t i = 0; i < gradients.size(); ++i)
+        jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+      static std::vector<FieldVector<double,1> > shapeValues;
+      localFE.localBasis().evaluateFunction(local, shapeValues);
+
+      static std::vector<FieldVector<double,1> > phaseShapeValues;
+      phaseLocalFE.localBasis().evaluateFunction(quadPoint.position(), phaseShapeValues);
+
+      double phase_atqp = 0.0;
+      for (std::size_t i = 0; i < phaseLocalFE.size(); ++i) {
+        auto row = phaseIndexSet.index(phaseLocalView.tree().localIndex(i));
+        phase_atqp += phase[row[0]] * phaseShapeValues[i];
+      }
+
+      for (size_t i=0; i<localFE.size(); i++) {
+        size_t row = localView.tree().localIndex(i);
+        for (size_t j=0; j<localFE.size(); j++ ) {
+          size_t col = localView.tree().localIndex(j);
+          elementMatrix[row][col] += (
+            phase_atqp * (gradients[i]*gradients[j])
+            + (1.0/(eps*eps*eps)) * (1.0-phase_atqp) * (shapeValues[i] * shapeValues[j]) ) * quadPoint.weight() * dx;
+        }
+
+        elementVector[row] += phase_atqp * shapeValues[i] * quadPoint.weight() * dx;
+      }
+    }
+  }
+}
+
+
+template <class Basis, class FineBasis, class Matrix, class Vector>
+void assembleSystem(const Basis& basis, const FineBasis& fine_basis, Vector const& phase, Matrix& matrix, Vector& rhs, double eps)
+{
+  using namespace Dune::Indices;
+
+  matrix.change_dim(basis.dimension(), basis.dimension());
+  matrix = 0;
+
+  rhs.change_dim(basis.dimension());
+  rhs = 0;
+
+  auto localView = basis.localView();
+  auto localIndexSet = basis.localIndexSet();
+
+  mtl::mat::inserter<Matrix, mtl::update_plus<double> > ins(matrix, 20);
+  for (const auto& element : elements(basis.gridView())) {
+    localView.bind(element);
+    localIndexSet.bind(localView);
+
+    mtl::dense2D<double> elementMatrix;
+    mtl::dense_vector<double> elementVector;
+    getLocalSystem(localView, fine_basis, phase, elementMatrix, elementVector, eps);