MultiBasis added and a system with for the diffuse-domain method

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
+Suggests: dune-uggrid dune-alugrid dune-foamgrid dune-functions dune-geometry

dune/multimesh/mmentity.hh

@@ -50,6 +50,21 @@ namespace Dune
       return LocalGeometry{LocalGeometryImp{(*this)[source_i], (*this)[target_i]}};
     }
 
+    /// Return coordinate `sourceLocal` in coordinates of target entity
+    static typename LocalGeometry::GlobalCoordinate local (
+        const HostGridEntity& source,
+        const HostGridEntity& target,
+        const typename LocalGeometry::LocalCoordinate& sourceLocal)
+    {
+      if (source.level() == target.level())
+        return sourceLocal;
+      else if (source.level() < target.level())
+        return localGeometry(target, source).local(sourceLocal);
+      else
+        return localGeometry(source, target).global(sourceLocal);
+    }
+
+
     /// \brief Return the entity seed which contains sufficient information
     /// to generate the entity again and uses as little memory as possible.
     /**

dune/multimesh/utility/CMakeLists.txt

 #install headers
 install(FILES
+  multibasis.hh
   structuredgridbuilder.hh
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/multimesh/utility)

dune/multimesh/utility/multibasis.hh

+#pragma once
+
+#include <tuple>
+#include <utility>
+
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/rangeutilities.hh>
+#include <dune/common/indices.hh>
+#include <dune/common/std/apply.hh>
+
+// requires the dune-functions module
+#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
+
+namespace Dune
+{
+  template <class MultiMesh, class... Bases>
+  struct MultiBasis
+  {
+    using MultiGridView = typename MultiMesh::LeafGridView;
+
+    /// Constructor. Stores a leafgridView of the MultiMesh and LocalView/LocalIndexSet of the bases
+    MultiBasis (MultiMesh const& mm, Bases&&... bases)
+      : multiGridView_(mm.leafGridView())
+      , bases_{std::forward<Bases>(bases)...}
+      , localViews_(Std::apply([](auto const&... b) { return std::make_tuple(b.localView()...); }, bases_))
+      , localIndexSets_(Std::apply([](auto const&... b) { return std::make_tuple(b.localIndexSet()...); }, bases_))
+    {}
+
+    /// Return the leafGridView of the MultiMesh
+    const MultiGridView& gridView () const { return multiGridView_; }
+
+    /// Return the I'th basis
+    template <std::size_t I>
+    auto const& basis (index_constant<I>) const { return std::get<I>(bases_); }
+
+    /// Return the LocalView assiciated with the I'th basis
+    template <std::size_t I>
+    auto const& localView (index_constant<I>) const { return std::get<I>(localViews_); }
+
+    /// Return the LocalIndexSet assiciated with the I'th basis
+    template <std::size_t I>
+    auto const& localIndexSet (index_constant<I>) const { return std::get<I>(localIndexSets_); }
+
+    /// Update all bases
+    void update ()
+    {
+      Hybrid::forEach(range(index_constant<sizeof...(Bases)>{}), [&,this](const auto i)
+      {
+        std::get<i.value>(bases_).update(multiGridView_.grid().leafGridView(i));
+      });
+    }
+
+    /// \brief Bind all localViews to the corresponding element and
+    /// bind the localIndexSets to these localViews
+    template <class MultiElement>
+    void bind (const MultiElement& multiElement)
+    {
+      Hybrid::forEach(range(index_constant<sizeof...(Bases)>{}), [&,this](const auto i)
+      {
+        std::get<i.value>(localViews_).bind(multiElement[i]);
+        std::get<i.value>(localIndexSets_).bind(std::get<i.value>(localViews_));
+      });
+    }
+
+  private:
+    MultiGridView multiGridView_;
+    std::tuple<Bases...> bases_;
+
+    std::tuple<typename Bases::LocalView...> localViews_;
+    std::tuple<typename Bases::LocalIndexSet...> localIndexSets_;
+  };
+
+  namespace Impl
+  {
+    template <class MM, class... Bases>
+    auto makeMultiBasisFinal(MM const& mm, Bases&&... bases)
+    {
+      return MultiBasis<MM, std::decay_t<Bases>...>{mm, std::forward<Bases>(bases)...};
+    }
+
+    template <class MM, class FactoryTuple, std::size_t... I>
+    auto makeMultiBasisImpl(MM const& mm, FactoryTuple&& factories, std::index_sequence<I...>)
+    {
+      using Functions::BasisBuilder::makeBasis;
+      return makeMultiBasisFinal(mm, makeBasis(mm[I].leafGridView(), std::get<I>(factories))...);
+    }
+  }
+
+  /// Construct a tuple of bases from its factories (factoryTags) and return a MultiBasis
+  template <class MM, class... Factories>
+  auto makeMultiBasis(MM const& mm, Factories&&... factories)
+  {
+    return Impl::makeMultiBasisImpl(mm, std::forward_as_tuple(factories...),
+      std::make_index_sequence<sizeof...(Factories)>{});
+  }
+
+} // end namespace Dune

src/CMakeLists.txt

@@ -9,3 +9,26 @@ target_link_dune_default_libraries("localrefinement")
 
 add_executable("uggrid" uggrid.cc)
 target_link_dune_default_libraries("uggrid")
+
+
+find_package(MTL PATHS /opt/development/mtl4)
+if (MTL_FOUND)
+  set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
+  set(MTL_COMPILE_DEFINITIONS "")
+  foreach(feature ${CXX_ELEVEN_FEATURE_LIST})
+    list(APPEND MTL_COMPILE_DEFINITIONS "MTL_WITH_${feature}")
+  endforeach()
+
+  if (HAVE_UMFPACK OR ENABLE_SUITESPARSE OR SuiteSparse_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)
+endif ()
\ No newline at end of file

src/macro.stand.2d

 DIM: 2
 DIM_OF_WORLD: 2
 
-number of elements: 4  
-number of vertices: 5  
+number of elements: 4
+number of vertices: 5
 
 element vertices:
-0 1 4 
-1 2 4  
-2 3 4 
-3 0 4 
+0 1 4
+1 2 4
+2 3 4
+3 0 4
 
 element boundaries:
-0 0 1 
+0 0 1
 0 0 2
-0 0 3 
-0 0 4 
+0 0 3
+0 0 4
 
 vertex coordinates:
- 0.0 0.0
- 1.0 0.0
- 1.0 1.0
- 0.0 1.0
- 0.5 0.5
+-1.5 -1.5
+ 1.5 -1.5
+ 1.5  1.5
+-1.5  1.5
+ 0.0  0.0
 
 element neighbours:
-1 3 -1 
-2 0 -1 
-3 1 -1 
+1 3 -1
+2 0 -1
+3 1 -1
 0 2 -1

src/phasefield.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/multimesh.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 <dune/multimesh/utility/multibasis.hh>
+#include "phasefield.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 HostGrid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
+  using Grid = MultiMesh<HostGrid>;
+
+  AlbertaReader<Grid> albertaReader;
+  GridFactory<Grid> factory(2);
+  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[0].globalRefine(6);
+  grid[1].globalRefine(4);
+
+  // refine grid[1] along phase-field interface
+  for (int i = 0; i < 10; ++i) {
+    std::cout << "refine " << i << "...\n";
+    for (auto const& elem : elements(grid[1].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,6-i)*eps;
+
+      if (refine)
+        grid[1].mark(1, elem);
+    }
+
+    grid[1].preAdapt();
+    grid[1].adapt();
+    grid[1].postAdapt();
+  }
+
+  using namespace Functions::BasisBuilder;
+  auto multiBasis = makeMultiBasis(grid, lagrange<1>(), lagrange<1>());
+
+  using VectorType = mtl::dense_vector<double>;
+  using MatrixType = mtl::compressed2D<double>;
+
+  // interpolate phase-field to fine grid
+  VectorType phase(multiBasis.basis(_1).dimension());
+  auto phaseWrapper = MTLVector(phase);
+  interpolate(multiBasis.basis(_1), phaseWrapper, [eps,d=signedDistFct](auto const& x)
+  {
+    return 0.5*(1.0 - std::tanh(3.0*d(x)/eps));
+  });
+  writeFile(phase, multiBasis.basis(_1), "phase");
+
+  // assemble a sequence of system for finer and finer grids
+  for (int i = 0; i < 6; ++i) {
+    grid[0].globalRefine(1);
+
+    VectorType rhs(multiBasis.basis(_0).dimension());
+    MatrixType matrix(multiBasis.basis(_0).dimension(), multiBasis.basis(_0).dimension());
+    assembleSystem(multiBasis, phase, matrix, rhs, eps);
+
+    VectorType u(multiBasis.basis(_0).dimension());
+    u = 0.0;
+    umfpack_solve(matrix, u, rhs);
+
+    writeFile(u, multiBasis.basis(_0), "u_" + std::to_string(i));
+  }
+}

src/phasefield.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 MultiBasis, class MultiElement>
+void getLocalSystem(const MultiBasis& multiBasis,
+                    const MultiElement& multiElement,
+                    mtl::dense_vector<double> const& phase,
+                    mtl::dense2D<double>& elementMatrix,
+                    mtl::dense_vector<double>& elementVector,
+                    double eps)
+{
+  using namespace Dune::Indices;
+
+  const int dim = MultiElement::dimension;
+  elementMatrix.change_dim(multiBasis.localView(_0).size(), multiBasis.localView(_0).size());
+  elementVector.change_dim(multiBasis.localView(_0).size());
+  set_to_zero(elementMatrix); // fills the entire matrix with zeroes
+  elementVector = 0;
+
+  // geometry of the basis element
+  auto geometry = multiElement[0].geometry();
+  const auto& localFE = multiBasis.localView(_0).tree().finiteElement();
+  const auto& phaseLocalFE = multiBasis.localView(_1).tree().finiteElement();
+
+  // the element to integrate over
+  auto const& quadElement = multiElement.max();
+  auto quadGeometry = quadElement.geometry();
+
+  const auto& quad = QuadratureRules<double, dim>::rule(quadGeometry.type(), 2*localFE.localBasis().order());
+  for (const auto& quadPoint : quad)
+  {
+    auto local = multiElement.local(quadElement, multiElement[0], quadPoint.position());
+    auto phaseLocal = multiElement.local(quadElement, multiElement[1], 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(phaseLocal, phaseShapeValues);
+
+    double phase_atqp = 0.0;
+    for (std::size_t i = 0; i < phaseLocalFE.size(); ++i) {
+      auto row = multiBasis.localIndexSet(_1).index(multiBasis.localView(_1).tree().localIndex(i));
+      phase_atqp += phase[row[0]] * phaseShapeValues[i];
+    }
+
+    for (size_t i=0; i<localFE.size(); i++) {
+      size_t row = multiBasis.localView(_0).tree().localIndex(i);
+      for (size_t j=0; j<localFE.size(); j++ ) {
+        size_t col = multiBasis.localView(_0).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] += shapeValues[i] * quadPoint.weight() * dx;
+    }
+  }
+}
+
+
+template <class MultiBasis, class Matrix, class Vector>
+void assembleSystem(MultiBasis& multiBasis, Vector const& phase, Matrix& matrix, Vector& rhs, double eps)
+{
+  using namespace Dune::Indices;
+
+  matrix.change_dim(multiBasis.basis(_0).dimension(), multiBasis.basis(_0).dimension());
+  matrix = 0;
+
+  rhs.change_dim(multiBasis.basis(_0).dimension());
+  rhs = 0;
+
+  multiBasis.update();
+
+  mtl::mat::inserter<Matrix, mtl::update_plus<double> > ins(matrix, 20);
+  for (const auto& multiElement : elements(multiBasis.gridView())) {
+    multiBasis.bind(multiElement);
+
+    mtl::dense2D<double> elementMatrix;
+    mtl::dense_vector<double> elementVector;
+    getLocalSystem(multiBasis, multiElement, phase, elementMatrix, elementVector, eps);
+
+    for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
+      auto row = multiBasis.localIndexSet(_0).index(i);
+      if (elementVector[i] != 0)
+        rhs[row[0]] += elementVector[i];
+
+      for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
+        auto col = multiBasis.localIndexSet(_0).index(j);
+        if (elementMatrix[i][j] != 0)
+          ins[row[0]][col[0]] += elementMatrix[i][j];
+      }
+    }
+  }
+}