Version 1.0

CMakeLists.txt

+# We require version CMake version 3.1 to prevent issues
+# with dune_enable_all_packages and older CMake versions.
+cmake_minimum_required(VERSION 3.1)
+project(dune-simon CXX)
+
+if(NOT (dune-common_DIR OR dune-common_ROOT OR
+      "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
+    string(REPLACE  ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR
+      ${PROJECT_BINARY_DIR})
+endif()
+
+#find dune-common and set the module path
+find_package(dune-common REQUIRED)
+list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules"
+  ${dune-common_MODULE_PATH})
+
+#include the dune macros
+include(DuneMacros)
+
+# start a dune project with information from dune.module
+dune_project()
+
+dune_enable_all_packages()
+
+add_subdirectory(src)
+add_subdirectory(dune)
+add_subdirectory(doc)
+add_subdirectory(cmake/modules)
+
+# finalize the dune project, e.g. generating config.h etc.
+finalize_dune_project(GENERATE_CONFIG_H_CMAKE)

README

+Preparing the Sources
+=========================
+
+Additional to the software mentioned in README you'll need the
+following programs installed on your system:
+
+  cmake >= 2.8.12
+
+Getting started
+---------------
+
+If these preliminaries are met, you should run
+
+  dunecontrol all
+
+which will find all installed dune modules as well as all dune modules
+(not installed) which sources reside in a subdirectory of the current
+directory. Note that if dune is not installed properly you will either
+have to add the directory where the dunecontrol script resides (probably
+./dune-common/bin) to your path or specify the relative path of the script.
+
+Most probably you'll have to provide additional information to dunecontrol
+(e. g. compilers, configure options) and/or make options.
+
+The most convenient way is to use options files in this case. The files
+define four variables:
+
+CMAKE_FLAGS      flags passed to cmake (during configure)
+
+An example options file might look like this:
+
+#use this options to configure and make if no other options are given
+CMAKE_FLAGS=" \
+-DCMAKE_CXX_COMPILER=g++-5 \
+-DCMAKE_CXX_FLAGS='-Wall -pedantic' \
+-DCMAKE_INSTALL_PREFIX=/install/path" #Force g++-5 and set compiler flags
+
+If you save this information into example.opts you can pass the opts file to
+dunecontrol via the --opts option, e. g.
+
+  dunecontrol --opts=example.opts all
+
+More info
+---------
+
+See
+
+     dunecontrol --help
+
+for further options.
+
+
+The full build system is described in the dune-common/doc/buildsystem (Git version) or under share/doc/dune-common/buildsystem if you installed DUNE!

cmake/modules/CMakeLists.txt

+set(modules "DuneSimonMacros.cmake")
+
+install(FILES ${modules} DESTINATION ${DUNE_INSTALL_MODULEDIR})

cmake/modules/DuneSimonMacros.cmake

+# File for module specific CMake tests.

config.h.cmake

+/* begin dune-simon
+   put the definitions for config.h specific to
+   your project here. Everything above will be
+   overwritten
+*/
+
+/* begin private */
+/* Name of package */
+#define PACKAGE "@DUNE_MOD_NAME@"
+
+/* Define to the address where bug reports for this package should be sent. */
+#define PACKAGE_BUGREPORT "@DUNE_MAINTAINER@"
+
+/* Define to the full name of this package. */
+#define PACKAGE_NAME "@DUNE_MOD_NAME@"
+
+/* Define to the full name and version of this package. */
+#define PACKAGE_STRING "@DUNE_MOD_NAME@ @DUNE_MOD_VERSION@"
+
+/* Define to the one symbol short name of this package. */
+#define PACKAGE_TARNAME "@DUNE_MOD_NAME@"
+
+/* Define to the home page for this package. */
+#define PACKAGE_URL "@DUNE_MOD_URL@"
+
+/* Define to the version of this package. */
+#define PACKAGE_VERSION "@DUNE_MOD_VERSION@"
+
+/* end private */
+
+/* Define to the version of dune-simon */
+#define DUNE_SIMON_VERSION "@DUNE_SIMON_VERSION@"
+
+/* Define to the major version of dune-simon */
+#define DUNE_SIMON_VERSION_MAJOR @DUNE_SIMON_VERSION_MAJOR@
+
+/* Define to the minor version of dune-simon */
+#define DUNE_SIMON_VERSION_MINOR @DUNE_SIMON_VERSION_MINOR@
+
+/* Define to the revision of dune-simon */
+#define DUNE_SIMON_VERSION_REVISION @DUNE_SIMON_VERSION_REVISION@
+
+/* end dune-simon
+   Everything below here will be overwritten
+*/

doc/CMakeLists.txt

+add_subdirectory("doxygen")

doc/doxygen/CMakeLists.txt

+# shortcut for creating the Doxyfile.in and Doxyfile
+add_doxygen_target()

doc/doxygen/Doxylocal

+# This file contains local changes to the doxygen configuration
+# please use '+=' to add files/directories to the lists
+
+# The INPUT tag can be used to specify the files and/or directories that contain
+# documented source files. You may enter file names like "myfile.cpp" or
+# directories like "/usr/src/myproject". Separate the files or directories
+# with spaces.
+
+INPUT                 += @top_srcdir@/dune/
+# see e.g. dune-grid for the examples of mainpage and modules
+# INPUT                 += @srcdir@/mainpage \
+#                          @srcdir@/modules
+
+# The EXCLUDE tag can be used to specify files and/or directories that should
+# be excluded from the INPUT source files. This way you can easily exclude a
+# subdirectory from a directory tree whose root is specified with the INPUT tag.
+
+# EXCLUDE               += @top_srcdir@/dune/simon/test
+
+# The EXAMPLE_PATH tag can be used to specify one or more files or
+# directories that contain example code fragments that are included (see
+# the \include command).
+
+# EXAMPLE_PATH          += @top_srcdir@/src
+
+# The IMAGE_PATH tag can be used to specify one or more files or
+# directories that contain image that are included in the documentation (see
+# the \image command).
+
+# IMAGE_PATH            += @top_srcdir@/dune/simon/pics

dune-simon.pc.in

+prefix=@prefix@
+exec_prefix=@exec_prefix@
+libdir=@libdir@
+includedir=@includedir@
+CXX=@CXX@
+CC=@CC@
+DEPENDENCIES=@REQUIRES@
+
+Name: @PACKAGE_NAME@
+Version: @VERSION@
+Description: dune-simon module
+URL: http://dune-project.org/
+Requires: dune-common dune-istl dune-typetree dune-geometry dune-localfunctions dune-uggrid dune-grid dune-functions
+Libs: -L${libdir}
+Cflags: -I${includedir}

dune.module

+################################
+# Dune module information file #
+################################
+
+# Name of the module
+Module: dune-simon
+Version: 1.0
+Maintainer: simon29996@mi.fu-berlin.de
+# Required build dependencies
+Depends: dune-common dune-istl dune-typetree dune-geometry dune-localfunctions dune-uggrid dune-grid dune-functions
+# Optional build dependencies
+#Suggests:

dune/CMakeLists.txt

+add_subdirectory(simon)

dune/simon/CMakeLists.txt

+#install headers
+install(FILES simon.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/simon)

dune/simon/simon.hh

+#ifndef DUNE_SIMON_HH
+#define DUNE_SIMON_HH
+
+// add your classes here
+
+#endif // DUNE_SIMON_HH

src/CMakeLists.txt

+add_executable("dune-simon" dune-simon.cc)
+target_link_dune_default_libraries("dune-simon")

src/dune-simon.cc

+#include <config.h>
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <map>
+#include <random>
+#include <cmath>
+
+#include <dune/common/timer.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/umfpack.hh>
+#include <dune/istl/cholmod.hh>
+#include <dune/istl/io.hh>
+#include <dune/istl/matrixmarket.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+using namespace Dune;
+
+
+void indexSet(const int& l, int& c, std::vector<int>& k0, std::vector<int>& k1)
+{
+  c = 0;
+
+  for(int i = 0; i <= l; i++){
+    int s;
+    if(i == 0){
+      s = 1;
+    }
+    else{
+      s = -l;
+    }
+
+    for(int j = s; j <= l; j++){
+      auto norm = std::sqrt(i*i + j*j);
+
+      if(norm <= l){
+        c += 1;
+        k0.push_back(i);
+        k1.push_back(j);
+      }
+    }
+  }
+}
+
+
+template <class VectorType>
+void generateSurfaceFunction(int& c, std::vector<int>& k0, std::vector<int>& k1, VectorType& xk,
+  VectorType& yk, const double& kappa, const double& sigma)
+{
+  double pi = std::acos(-1.0);
+
+  xk.resize(c);
+  yk.resize(c);
+
+  std::random_device rd{};
+  std::mt19937 gen{rd()};
+  std::normal_distribution<> d;
+
+  for(int i = 0; i < c; i++){
+    auto normsq = k0[i]*k0[i] + k1[i]*k1[i];
+    double factor = std::sqrt(2*normsq*(4*pi*pi*kappa*normsq + sigma));
+    xk[i] = d(gen)/factor;
+    yk[i] = d(gen)/factor;
+  }
+}
+
+
+template <class CoordinateType, class VectorType, class MatrixType>
+void getFundamentalMatrix(CoordinateType& x, MatrixType& fundamentalmatrix, int& c,
+  std::vector<int>& k0, std::vector<int>& k1, VectorType& xk, VectorType& yk, double& pi)
+{
+  double sp;
+  double factor;
+  std::vector<double> gradient;
+  gradient.resize(2);
+  gradient[0] = 0;
+  gradient[1] = 0;
+
+  for(int i=0; i < c; i++){
+    sp = 2*pi*(x[0]*k0[i] + x[1]*k1[i]);
+    factor = std::cos(sp)*xk[i] + std::sin(sp)*yk[i];
+
+    gradient[0] += k0[i] * factor;
+    gradient[1] += k1[i] * factor;
+  }
+
+  fundamentalmatrix[0][0] = 1 + 4*gradient[0]*gradient[0];
+  fundamentalmatrix[0][1] = 4*gradient[0]*gradient[1];
+  fundamentalmatrix[1][0] = 4*gradient[1]*gradient[0];
+  fundamentalmatrix[1][1] = 1 + 4*gradient[1]*gradient[1];
+}
+
+
+template <class LocalView, class VectorType, class MatrixType>
+void getLocalMatrices(const LocalView& localview, MatrixType& localstiffnessmatrix,
+  MatrixType& localmassmatrix, int& c, std::vector<int>& k0, std::vector<int>& k1,
+  VectorType& xk, VectorType& yk, double& pi)
+{
+  // Get the grid element from the local FE basis view
+  typedef typename LocalView::Element Element;
+  const Element& element = localview.element();
+
+  const int dim = Element::dimension;
+  auto geometry = element.geometry();
+
+  // Get set of shape functions for this element
+  const auto& localFiniteElement = localview.tree().finiteElement();
+
+  // Set all matrix entries to zero
+  localstiffnessmatrix.setSize(localFiniteElement.localBasis().size(), localFiniteElement.localBasis().size());
+  localstiffnessmatrix = 0;
+  localmassmatrix.setSize(localFiniteElement.localBasis().size(), localFiniteElement.localBasis().size());
+  localmassmatrix = 0;
+
+  // Get a quadrature rule
+  const int quadOrder = dim*localFiniteElement.localBasis().order();
+  const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), quadOrder);
+
+  // Loop over all quadrature points
+  for (size_t pt=0; pt < quad.size(); pt++) {
+
+    // Position of the current quadrature point in the reference element
+    const FieldVector<double,dim>& quadPos = quad[pt].position();
+    // std::cout << quadPos << std::endl;
+
+    // Position of the current quadrature point in the actual element
+    const auto globalQuadPos = geometry.global(quadPos);
+
+    // The transposed inverse Jacobian of the map from the reference element to the element
+    const auto& jacobian = geometry.jacobianInverseTransposed(quadPos);
+
+    // The multiplicative factor in the integral transformation formula
+    const double integrationElement = geometry.integrationElement(quadPos);
+
+    // The gradients of the shape functions on the reference element
+    std::vector<FieldMatrix<double,1,dim>> referenceGradients;
+    localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+    // Compute the shape function gradients on the real element
+    std::vector<FieldVector<double,dim>> gradients(referenceGradients.size());
+    for (size_t i=0; i<gradients.size(); i++){
+      jacobian.mv(referenceGradients[i][0], gradients[i]);
+    }
+
+    // The values of the shape functions on the reference element
+    std::vector<FieldVector<double,1>> values;
+    localFiniteElement.localBasis().evaluateFunction(quadPos, values);
+
+    FieldMatrix<double,dim,dim> fundamentalMatrix;
+    getFundamentalMatrix(globalQuadPos, fundamentalMatrix, c, k0, k1, xk, yk, pi);
+    auto transformationFactor = std::sqrt(fundamentalMatrix.determinant());
+    auto transformationMatrix = fundamentalMatrix;
+    transformationMatrix.invert();
+
+    // Compute the actual matrix entries
+    FieldVector<double,dim> x;
+
+    for (size_t i=0; i<localstiffnessmatrix.N(); i++){
+
+      transformationMatrix.mv(gradients[i],x);
+
+      for (size_t j=0; j<localstiffnessmatrix.M(); j++ ){
+        localstiffnessmatrix[i][j] += transformationFactor * ( x*gradients[j] ) * quad[pt].weight() * integrationElement;
+        localmassmatrix[i][j] += transformationFactor * (values[i]*values[j]) * quad[pt].weight() * integrationElement;
+      }
+    }
+  }
+
+}
+
+
+template <class IndexType>
+bool isOnBoundary(IndexType& index, const int& m, const int& feorder)
+{
+  // Technicality to avoid comparing int with unsigned int
+  int Index = index;
+
+  if(Index <= m*(m+2)){
+    if(Index >= (m+1)*m){
+      return true;
+    }
+    else if(Index%(m+1) == m){
+      return true;
+    }
+  }
+  else if(Index <= m*((m+1)*feorder + 1)){
+    if(Index >= m*(m*feorder + 2) + 1){
+      return true;
+    }
+  }
+  else if(Index >= (m+1)*(m*feorder + 1) && Index <= m*((m+1)*(2*feorder - 1) + 1)){
+    if((Index - (m+1)*(m+1))%((m+1)*(feorder-1)) >= m*(feorder - 1)){
+      return true;
+    }
+  }
+  return false;
+}
+
+
+template <class IndexType>
+void indexCorrection(IndexType& index, const int& m, const int& feorder)
+{
+  // Technicality to avoid comparing int with unsigned int
+  int Index = index;
+
+  if(Index <= m*(m+2)){
+    if(Index == m*(m+2)){
+      Index = 0;
+    }
+    else if(Index >= (m+1)*m){
+      Index -= (m+1)*m;
+    }
+    else if(Index%(m+1) == m){
+      Index -= m;
+    }
+  }
+  else if(Index <= m*((m+1)*feorder + 1)){
+    if(Index >= m*(m*feorder + 2) + 1){
+      Index -= m*m*(feorder-1);
+    }
+  }
+  else if(Index >= (m+1)*(m*feorder + 1) && Index <= m*((m+1)*(2*feorder - 1) + 1)){
+    if((Index - (m+1)*(m+1))%((m+1)*(feorder-1)) >= m*(feorder - 1)){
+      Index -= m*(feorder-1);
+    }
+  }
+
+  index = Index;
+}
+
+
+// Get the occupation pattern of the stiffness/mass matrix
+template <class FEBasis>
+void getOccupationPattern(const FEBasis& febasis, MatrixIndexSet& matrixindexset, const int& m,
+  const int& feorder)
+{
+  // Total number of grid vertices
+  const int n = febasis.size();
+
+  matrixindexset.resize(n, n);
+
+  // A view on the FE basis on a single element
+  auto localView = febasis.localView();
+
+  // Loop over all leaf elements
+  for(const auto& e : elements(febasis.gridView()))
+  {
+    // Bind the local FE basis view to the current element
+    localView.bind(e);
+
+    // There is a matrix entry a_ij if the i-th and j-th vertex are connected in the grid
+    for (size_t i=0; i<localView.tree().size(); i++) {
+
+      auto iIdx = localView.index(i);
+      indexCorrection(iIdx[0], m, feorder);
+
+      for (size_t j=0; j<localView.tree().size(); j++) {
+
+        auto jIdx = localView.index(j);
+        indexCorrection(jIdx[0], m, feorder);
+
+        // Add a nonzero entry to the matrix
+        matrixindexset.add(iIdx, jIdx);
+
+      }
+
+    }
+
+  }
+
+  for(int i = 0; i < n; i++){
+    if(isOnBoundary(i, m, feorder)){
+      matrixindexset.add(i,i);
+    }
+  }
+
+}
+
+
+//Assemble the Laplace stiffness and mass matrix on the given grid view
+template <class FEBasis, class VectorType, class MatrixType>
+void assembleMatrices(const FEBasis& febasis, MatrixType& stiffnessmatrix, MatrixType& massmatrix,
+  const int& m, const int& feorder, int& c, std::vector<int>& k0, std::vector<int>& k1, VectorType& xk,
+  VectorType& yk, double& pi)
+{
+  const int n = febasis.size();
+  stiffnessmatrix.setSize(n,n);
+  massmatrix.setSize(n,n);
+
+  // Get the grid view from the finite element basis
+  typedef typename FEBasis::GridView GridView;
+  GridView gridView = febasis.gridView();
+
+  // MatrixIndexSets store the occupation pattern of a sparse matrix.
+  // They are not particularly efficient, but simple to use.
+  MatrixIndexSet occupationPattern;
+  getOccupationPattern(febasis, occupationPattern, m, feorder);
+
+  // ... and give it the occupation pattern we want.
+  occupationPattern.exportIdx(stiffnessmatrix);
+  occupationPattern.exportIdx(massmatrix);
+
+  // Set all entries to zero
+  stiffnessmatrix = 0;
+  massmatrix = 0;
+
+  // A view on the FE basis on a single element
+  auto localView = febasis.localView();
+
+  // A loop over all elements of the grid
+  for(const auto& e : elements(gridView))
+  {
+
+    // Bind the local FE basis view to the current element
+    localView.bind(e);
+
+    // Now let's get the element stiffness matrix
+    // A dense matrix is used for the element stiffness matrix
+    Matrix<FieldMatrix<double,1,1>> localStiffnessMatrix;
+    // getLocalStiffnessMatrix(localView, localstiffnessmatrix, c, k0, k1, xk, yk);
+    Matrix<FieldMatrix<double,1,1>> localMassMatrix;
+    // getLocalMassMatrix(localView, localmassmatrix, c, k0, k1, xk, yk);
+    getLocalMatrices(localView, localStiffnessMatrix, localMassMatrix, c, k0, k1, xk, yk, pi);
+
+    // Add element stiffness matrix onto the global stiffness matrix
+    for (size_t i=0; i<localStiffnessMatrix.N(); i++) {
+
+      // The global index of the i-th local degree of freedom of the element 'e'
+      auto row = localView.index(i);
+      indexCorrection(row[0], m, feorder);
+
+      for (size_t j=0; j<localStiffnessMatrix.M(); j++ ) {
+
+        // The global index of the j-th local degree of freedom of the element 'e'
+        auto col = localView.index(j);
+        indexCorrection(col[0], m, feorder);
+
+        stiffnessmatrix[row][col] += localStiffnessMatrix[i][j];
+        massmatrix[row][col] += localMassMatrix[i][j];
+
+      }
+
+    }
+
+  }
+
+  for(int i = 0; i < n; i++){
+    if(isOnBoundary(i, m, feorder)){
+      stiffnessmatrix[i][i] = 1;
+    }
+  }
+
+}
+
+
+template <class VectorType>
+void copyToBoundary(VectorType& vector, const int& n, const int& m, const int& feorder)
+{
+  for(int i = 0; i < n; i++){
+    if(isOnBoundary(i, m, feorder)){
+      auto j = i;
+      indexCorrection(j, m, feorder);
+      vector[i] = vector[j];
+    }
+  }
+}
+
+
+template <class FEBasis, class GridView, class VectorType>
+void interpolateSurfaceFunction(const FEBasis& febasis, GridView& gridview, VectorType& surface,
+  std::vector<char>& nonboundarynodes, const int& m, const int& feorder, int& c, std::vector<int>& k0,
+  std::vector<int>& k1, VectorType& xk, VectorType& yk, double& pi)
+{
+  const int n = febasis.size();
+
+  auto graphFunction  = [c, k0, k1, xk, yk, pi](const auto x) {
+    double value = 0;
+    double sp;
+
+    for(int i=0; i < c; i++){
+      sp = 2*pi*(x[0]*k0[i] + x[1]*k1[i]);
+      value += std::sin(sp)*xk[i] - std::cos(sp)*yk[i];
+    }
+    value /= pi;
+    return value;
+  };
+
+  surface.resize(n);
+
+  interpolate(febasis, surface, graphFunction, nonboundarynodes);
+  copyToBoundary(surface, n, m, feorder);
+}
+
+
+template <class FEBasis, class Solver, class VectorType, class MatrixType>
+void implicitEuler(const FEBasis& febasis, Solver& solver, VectorType& iter,
+  MatrixType& massmatrix, const int& m)
+{
+  const int n = febasis.size();
+
+  VectorType rhs;
+  rhs.resize(n);
+  rhs = 0;
+
+  massmatrix.mv(iter,rhs);
+
+  iter = 0;
+
+  InverseOperatorResult statistics;
+  solver.apply(iter, rhs, statistics);
+}
+
+
+int main (int argc, char *argv[]) try
+{
+  MPIHelper::instance(argc, argv);
+
+
+  // Initialize parameters
+  const int m = 400;
+  const int feOrder = 4;
+  const int l = 100;
+  const double kappa = 0;
+  const double sigma = 1;
+
+  const int refinementNumber = 2;
+
+  const double stepSize = 0.0005;
+  const int steps = 100;
+  const int skipExport = 100;
+
+
+  // Initialize options
+  std::string answer1;
+  std::cout << "Do you want to generate a new surface function? (y/n) ";
+  std::cin >> answer1;
+
+  std::string answer2;
+  if(answer1 == "y"){
+    std::cout << "Do you want to save the current surface function? (y/n) ";
+    std::cin >> answer2;
+  }
+
+  std::string answer3;
+  std::cout << "Do you want to generate new matrices? (y/n) ";
+  std::cin >> answer3;
+
+  std::string answer4;
+  if(answer3 == "y"){
+    std::cout << "Do you want to save the current matrices? (y/n) ";
+    std::cin >> answer4;
+  }
+
+  std::string answer5;
+  std::cout << "Do you want to export the surface? (y/n) ";
+  std::cin >> answer5;
+
+
+  // Initialize other stuff
+  std::string name1;
+  std::string name2;
+  name1 = "_" + std::to_string(l) + "_" + std::to_string(kappa) + "_" + std::to_string(sigma);
+  name2 = name1 + "_" + std::to_string(m) + "_" + std::to_string(feOrder);
+
+  typedef BlockVector<FieldVector<double,1>> VectorType;
+  typedef BCRSMatrix<FieldMatrix<double,1,1>> MatrixType;
+
+  double pi = std::acos(-1.0);
+
+  Timer timer;
+
+
+  // Setting up the random surface
+  int c;
+  std::vector<int> k0;
+  std::vector<int> k1;
+
+  indexSet(l,c,k0,k1);
+
+  std::cout << "Number of Fourier Nodes is " << 2*c << std::endl;
+
+  VectorType Xk;
+  VectorType Yk;
+
+  if(answer1 == "y"){
+    generateSurfaceFunction(c, k0, k1, Xk, Yk, kappa, sigma);
+
+    if(answer2 == "y"){
+      storeMatrixMarket(Xk, "X" + name1);
+      storeMatrixMarket(Yk, "Y" + name1);
+    }
+  }
+  else{
+    loadMatrixMarket(Xk, "X" + name1);
+    loadMatrixMarket(Yk, "Y" + name1);
+  }
+
+
+  // Generation of the grid
+  Dune::YaspGrid<2> grid({1,1},{m,m});
+
+  auto gridView = grid.leafGridView();
+
+  using GridView = decltype(gridView);
+
+
+  // Assembling of the matrices
+  typedef Functions::LagrangeBasis<GridView,feOrder> FEBasis;
+  FEBasis feBasis(gridView);
+
+  auto localView = feBasis.localView();
+
+  const int n = feBasis.size();
+  std::cout << "Number of DOFs is " << n << std::endl;
+
+  MatrixType stiffnessMatrix;
+  MatrixType massMatrix;
+
+  if(answer3 == "y"){
+    timer.reset();
+    assembleMatrices(feBasis, stiffnessMatrix, massMatrix, m, feOrder, c, k0, k1, Xk, Yk, pi);
+    std::cout << "Assembling the matrices took " << timer.elapsed() << "s" << std::endl;
+
+    if(answer4 == "y"){
+      storeMatrixMarket(stiffnessMatrix,"stiffnessmatrix" + name2);
+      storeMatrixMarket(massMatrix,"massmatrix" + name2);
+    }
+  }
+  else{
+    loadMatrixMarket(stiffnessMatrix,"stiffnessmatrix" + name2);
+    loadMatrixMarket(massMatrix,"massmatrix" + name2);
+  }
+
+  auto& operatorMatrix = stiffnessMatrix;
+  operatorMatrix *= stepSize;
+  operatorMatrix += massMatrix;
+
+
+  // Save the nonboundary nodes
+  std::vector<char> nonBoundaryNodes;
+  nonBoundaryNodes.resize(n);
+  for(int i = 0; i < n; i++){
+    if(isOnBoundary(i, m, feOrder) == false){
+      nonBoundaryNodes[i] = true;
+    }
+  }
+
+
+  // Generating the surface data
+  VectorType surface;
+  if(answer5 == "y"){
+    timer.reset();
+    interpolateSurfaceFunction(feBasis, gridView, surface, nonBoundaryNodes, m, feOrder, c, k0, k1, Xk, Yk, pi);
+    std::cout << "Interpolating the surface function took " << timer.elapsed() << "s" << std::endl;
+  }
+
+
+  // Initial value handling
+  VectorType iter;
+  iter.resize(n);
+
+  auto initialValueFunction  = [pi](const auto x) { return 10000/pi*std::exp(-10000*((x[0]-0.5)*(x[0]-0.5) + (x[1]-0.5)*(x[1]-0.5))); };
+  interpolate(feBasis, iter, initialValueFunction, nonBoundaryNodes);
+  copyToBoundary(iter, n, m, feOrder);
+
+
+  // Generating the output files
+  auto iterFunction = Functions::makeDiscreteGlobalBasisFunction<double>(feBasis, iter);
+  auto surfaceFunction = Functions::makeDiscreteGlobalBasisFunction<double>(feBasis, surface);
+
+  auto vtkWriter = std::make_shared<SubsamplingVTKWriter<GridView>>(gridView, refinementLevels(refinementNumber));
+  VTKSequenceWriter<GridView> vtkSequenceWriter(vtkWriter, "solution");
+  vtkWriter->addVertexData(iterFunction, VTK::FieldInfo("solution", VTK::FieldInfo::Type::scalar, 1));
+  if(answer5 == "y"){
+    vtkWriter->addVertexData(surfaceFunction, VTK::FieldInfo("surface", VTK::FieldInfo::Type::scalar, 1));
+  }
+  vtkSequenceWriter.write(0.0);
+
+
+  // Time integration
+  double time = 0;
+
+  // Setting up the solver
+  timer.reset();
+  // UMFPack<MatrixType> directSolver(operatorMatrix, false);
+  Cholmod<VectorType> directSolver;
+  directSolver.setMatrix(operatorMatrix);
+  std::cout << "Solving the linear system took " << timer.elapsed() << "s" << std::endl;
+
+  for (int i = 1; i <= steps; i++){
+
+    time += stepSize;
+
+    timer.reset();
+    implicitEuler(feBasis, directSolver, iter, massMatrix, m);
+    std::cout << timer.elapsed() << "s" << std::endl;
+
+    if (i % skipExport == 0){
+      copyToBoundary(iter, n, m, feOrder);
+      vtkSequenceWriter.write(time);
+    }
+
+  }
+
+
+  return 0;
+}
+
+
+// Error handling
+catch (Exception& e) {
+  std::cout << e.what() << std::endl;
+}