dune/tectonic/time-stepping/state/ageinglawstateupdater.hh

+#ifndef SRC_TIME_STEPPING_STATE_AGEINGLAWSTATEUPDATER_HH
+#define SRC_TIME_STEPPING_STATE_AGEINGLAWSTATEUPDATER_HH
+
+#include <dune/common/bitsetvector.hh>
+
+#include "stateupdater.hh"
+
+template <class ScalarVector, class Vector>
+class AgeingLawStateUpdater : public LocalStateUpdater<ScalarVector, Vector> {
+private:
+    using BitVector = Dune::BitSetVector<1>;
+
+public:
+    AgeingLawStateUpdater(
+            const ScalarVector& alpha_initial,
+            const BitVector& nodes,
+            const double L,
+            const double V0);
+
+  void nextTimeStep() override;
+  void setup(double tau) override;
+  void solve(const Vector& velocity_field) override;
+  void extractAlpha(ScalarVector&) override;
+
+  auto clone() const -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> override;
+
+private:
+  std::vector<int> localToGlobal_;
+
+  ScalarVector alpha_o_;
+  ScalarVector alpha_;
+  const BitVector& nodes_;
+  const double L_;
+  const double V0_;
+  double tau_;
+};
+#endif

dune/tectonic/time-stepping/state/explicit.aux

+\relax 

dune/tectonic/time-stepping/state/explicit.log

+This is pdfTeX, Version 3.14159265-2.6-1.40.17 (TeX Live 2016/Debian) (preloaded format=pdflatex 2018.9.22)  2 SEP 2019 17:19
+entering extended mode
+ restricted \write18 enabled.
+ %&-line parsing enabled.
+**explicit.tex
+(./explicit.tex
+LaTeX2e <2017/01/01> patch level 3
+Babel <3.9r> and hyphenation patterns for 83 language(s) loaded.
+(/usr/share/texlive/texmf-dist/tex/latex/standalone/standalone.cls
+Document Class: standalone 2015/07/15 v1.2 Class to compile TeX sub-files stand
+alone
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ifluatex.sty
+Package: ifluatex 2016/05/16 v1.4 Provides the ifluatex switch (HO)
+Package ifluatex Info: LuaTeX not detected.
+)
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ifpdf.sty
+Package: ifpdf 2016/05/14 v3.1 Provides the ifpdf switch
+)
+(/usr/share/texlive/texmf-dist/tex/generic/ifxetex/ifxetex.sty
+Package: ifxetex 2010/09/12 v0.6 Provides ifxetex conditional
+)
+(/usr/share/texlive/texmf-dist/tex/latex/xkeyval/xkeyval.sty
+Package: xkeyval 2014/12/03 v2.7a package option processing (HA)
+
+(/usr/share/texlive/texmf-dist/tex/generic/xkeyval/xkeyval.tex
+(/usr/share/texlive/texmf-dist/tex/generic/xkeyval/xkvutils.tex
+\XKV@toks=\toks14
+\XKV@tempa@toks=\toks15
+
+(/usr/share/texlive/texmf-dist/tex/generic/xkeyval/keyval.tex))
+\XKV@depth=\count79
+File: xkeyval.tex 2014/12/03 v2.7a key=value parser (HA)
+))
+\sa@internal=\count80
+\c@sapage=\count81
+
+(/usr/share/texlive/texmf-dist/tex/latex/standalone/standalone.cfg
+File: standalone.cfg 2015/07/15 v1.2 Default configuration file for 'standalone
+' class
+)
+(/usr/share/texlive/texmf-dist/tex/latex/base/article.cls
+Document Class: article 2014/09/29 v1.4h Standard LaTeX document class
+(/usr/share/texlive/texmf-dist/tex/latex/base/size10.clo
+File: size10.clo 2014/09/29 v1.4h Standard LaTeX file (size option)
+)
+\c@part=\count82
+\c@section=\count83
+\c@subsection=\count84
+\c@subsubsection=\count85
+\c@paragraph=\count86
+\c@subparagraph=\count87
+\c@figure=\count88
+\c@table=\count89
+\abovecaptionskip=\skip41
+\belowcaptionskip=\skip42
+\bibindent=\dimen102
+)
+(/usr/share/texmf/tex/latex/preview/preview.sty
+Package: preview 2010/02/14 11.90 (AUCTeX/preview-latex)
+
+(/usr/share/texmf/tex/latex/preview/prtightpage.def
+\PreviewBorder=\dimen103
+)
+\pr@snippet=\count90
+\pr@box=\box26
+\pr@output=\toks16
+))
+(/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsmath.sty
+Package: amsmath 2016/11/05 v2.16a AMS math features
+\@mathmargin=\skip43
+
+For additional information on amsmath, use the `?' option.
+(/usr/share/texlive/texmf-dist/tex/latex/amsmath/amstext.sty
+Package: amstext 2000/06/29 v2.01 AMS text
+
+(/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsgen.sty
+File: amsgen.sty 1999/11/30 v2.0 generic functions
+\@emptytoks=\toks17
+\ex@=\dimen104
+))
+(/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsbsy.sty
+Package: amsbsy 1999/11/29 v1.2d Bold Symbols
+\pmbraise@=\dimen105
+)
+(/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsopn.sty
+Package: amsopn 2016/03/08 v2.02 operator names
+)
+\inf@bad=\count91
+LaTeX Info: Redefining \frac on input line 213.
+\uproot@=\count92
+\leftroot@=\count93
+LaTeX Info: Redefining \overline on input line 375.
+\classnum@=\count94
+\DOTSCASE@=\count95
+LaTeX Info: Redefining \ldots on input line 472.
+LaTeX Info: Redefining \dots on input line 475.
+LaTeX Info: Redefining \cdots on input line 596.
+\Mathstrutbox@=\box27
+\strutbox@=\box28
+\big@size=\dimen106
+LaTeX Font Info:    Redeclaring font encoding OML on input line 712.
+LaTeX Font Info:    Redeclaring font encoding OMS on input line 713.
+\macc@depth=\count96
+\c@MaxMatrixCols=\count97
+\dotsspace@=\muskip10
+\c@parentequation=\count98
+\dspbrk@lvl=\count99
+\tag@help=\toks18
+\row@=\count100
+\column@=\count101
+\maxfields@=\count102
+\andhelp@=\toks19
+\eqnshift@=\dimen107
+\alignsep@=\dimen108
+\tagshift@=\dimen109
+\tagwidth@=\dimen110
+\totwidth@=\dimen111
+\lineht@=\dimen112
+\@envbody=\toks20
+\multlinegap=\skip44
+\multlinetaggap=\skip45
+\mathdisplay@stack=\toks21
+LaTeX Info: Redefining \[ on input line 2817.
+LaTeX Info: Redefining \] on input line 2818.
+)
+(/usr/share/texlive/texmf-dist/tex/latex/mathtools/mathtools.sty
+Package: mathtools 2015/11/12 v1.18 mathematical typesetting tools
+
+(/usr/share/texlive/texmf-dist/tex/latex/tools/calc.sty
+Package: calc 2014/10/28 v4.3 Infix arithmetic (KKT,FJ)
+\calc@Acount=\count103
+\calc@Bcount=\count104
+\calc@Adimen=\dimen113
+\calc@Bdimen=\dimen114
+\calc@Askip=\skip46
+\calc@Bskip=\skip47
+LaTeX Info: Redefining \setlength on input line 80.
+LaTeX Info: Redefining \addtolength on input line 81.
+\calc@Ccount=\count105
+\calc@Cskip=\skip48
+)
+(/usr/share/texlive/texmf-dist/tex/latex/mathtools/mhsetup.sty
+Package: mhsetup 2010/01/21 v1.2a programming setup (MH)
+)
+LaTeX Info: Thecontrolsequence`\('isalreadyrobust on input line 129.
+LaTeX Info: Thecontrolsequence`\)'isalreadyrobust on input line 129.
+LaTeX Info: Thecontrolsequence`\['isalreadyrobust on input line 129.
+LaTeX Info: Thecontrolsequence`\]'isalreadyrobust on input line 129.
+\g_MT_multlinerow_int=\count106
+\l_MT_multwidth_dim=\dimen115
+\origjot=\skip49
+\l_MT_shortvdotswithinadjustabove_dim=\dimen116
+\l_MT_shortvdotswithinadjustbelow_dim=\dimen117
+\l_MT_above_intertext_sep=\dimen118
+\l_MT_below_intertext_sep=\dimen119
+\l_MT_above_shortintertext_sep=\dimen120
+\l_MT_below_shortintertext_sep=\dimen121
+)
+No file explicit.aux.
+\openout1 = `explicit.aux'.
+
+LaTeX Font Info:    Checking defaults for OML/cmm/m/it on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+LaTeX Font Info:    Checking defaults for T1/cmr/m/n on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+LaTeX Font Info:    Checking defaults for OT1/cmr/m/n on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+LaTeX Font Info:    Checking defaults for OMS/cmsy/m/n on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+LaTeX Font Info:    Checking defaults for OMX/cmex/m/n on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+LaTeX Font Info:    Checking defaults for U/cmr/m/n on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+Preview: Fontsize 10pt
+Preview: PDFoutput 1
+(/usr/share/texlive/texmf-dist/tex/latex/graphics/graphicx.sty
+Package: graphicx 2014/10/28 v1.0g Enhanced LaTeX Graphics (DPC,SPQR)
+
+(/usr/share/texlive/texmf-dist/tex/latex/graphics/graphics.sty
+Package: graphics 2016/10/09 v1.0u Standard LaTeX Graphics (DPC,SPQR)
+
+(/usr/share/texlive/texmf-dist/tex/latex/graphics/trig.sty
+Package: trig 2016/01/03 v1.10 sin cos tan (DPC)
+)
+(/usr/share/texlive/texmf-dist/tex/latex/graphics-cfg/graphics.cfg
+File: graphics.cfg 2016/06/04 v1.11 sample graphics configuration
+)
+Package graphics Info: Driver file: pdftex.def on input line 99.
+
+(/usr/share/texlive/texmf-dist/tex/latex/graphics-def/pdftex.def
+File: pdftex.def 2017/01/12 v0.06k Graphics/color for pdfTeX
+
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/infwarerr.sty
+Package: infwarerr 2016/05/16 v1.4 Providing info/warning/error messages (HO)
+)
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ltxcmds.sty
+Package: ltxcmds 2016/05/16 v1.23 LaTeX kernel commands for general use (HO)
+)
+\Gread@gobject=\count107
+
+(/usr/share/texlive/texmf-dist/tex/context/base/mkii/supp-pdf.mkii
+[Loading MPS to PDF converter (version 2006.09.02).]
+\scratchcounter=\count108
+\scratchdimen=\dimen122
+\scratchbox=\box29
+\nofMPsegments=\count109
+\nofMParguments=\count110
+\everyMPshowfont=\toks22
+\MPscratchCnt=\count111
+\MPscratchDim=\dimen123
+\MPnumerator=\count112
+\makeMPintoPDFobject=\count113
+\everyMPtoPDFconversion=\toks23
+))) (/usr/share/texlive/texmf-dist/tex/generic/oberdiek/pdftexcmds.sty
+Package: pdftexcmds 2016/05/21 v0.22 Utility functions of pdfTeX for LuaTeX (HO
+)
+Package pdftexcmds Info: LuaTeX not detected.
+Package pdftexcmds Info: \pdf@primitive is available.
+Package pdftexcmds Info: \pdf@ifprimitive is available.
+Package pdftexcmds Info: \pdfdraftmode found.
+)
+(/usr/share/texlive/texmf-dist/tex/latex/oberdiek/epstopdf-base.sty
+Package: epstopdf-base 2016/05/15 v2.6 Base part for package epstopdf
+
+(/usr/share/texlive/texmf-dist/tex/latex/oberdiek/grfext.sty
+Package: grfext 2016/05/16 v1.2 Manage graphics extensions (HO)
+
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/kvdefinekeys.sty
+Package: kvdefinekeys 2016/05/16 v1.4 Define keys (HO)
+))
+(/usr/share/texlive/texmf-dist/tex/latex/oberdiek/kvoptions.sty
+Package: kvoptions 2016/05/16 v3.12 Key value format for package options (HO)
+
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/kvsetkeys.sty
+Package: kvsetkeys 2016/05/16 v1.17 Key value parser (HO)
+
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/etexcmds.sty
+Package: etexcmds 2016/05/16 v1.6 Avoid name clashes with e-TeX commands (HO)
+Package etexcmds Info: Could not find \expanded.
+(etexcmds)             That can mean that you are not using pdfTeX 1.50 or
+(etexcmds)             that some package has redefined \expanded.
+(etexcmds)             In the latter case, load this package earlier.
+)))
+Package epstopdf-base Info: Redefining graphics rule for `.eps' on input line 4
+38.
+Package grfext Info: Graphics extension search list:
+(grfext)             [.png,.pdf,.jpg,.mps,.jpeg,.jbig2,.jb2,.PNG,.PDF,.JPG,.JPE
+G,.JBIG2,.JB2,.eps]
+(grfext)             \AppendGraphicsExtensions on input line 456.
+
+(/usr/share/texlive/texmf-dist/tex/latex/latexconfig/epstopdf-sys.cfg
+File: epstopdf-sys.cfg 2010/07/13 v1.3 Configuration of (r)epstopdf for TeX Liv
+e
+))
+\Gin@req@height=\dimen124
+\Gin@req@width=\dimen125
+)
+Preview: Tightpage 0 0 0 0
+[1{/var/lib/texmf/fonts/map/pdftex/updmap/pdftex.map}] (./explicit.aux) ) 
+Here is how much of TeX's memory you used:
+ 3166 strings out of 493013
+ 46673 string characters out of 6135682
+ 112782 words of memory out of 5000000
+ 6679 multiletter control sequences out of 15000+600000
+ 4094 words of font info for 16 fonts, out of 8000000 for 9000
+ 1141 hyphenation exceptions out of 8191
+ 53i,12n,56p,255b,147s stack positions out of 5000i,500n,10000p,200000b,80000s
+</usr/
+share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmex10.pfb></usr/share/
+texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmmi10.pfb></usr/share/texliv
+e/texmf-dist/fonts/type1/public/amsfonts/cm/cmmi5.pfb></usr/share/texlive/texmf
+-dist/fonts/type1/public/amsfonts/cm/cmmi7.pfb></usr/share/texlive/texmf-dist/f
+onts/type1/public/amsfonts/cm/cmr10.pfb></usr/share/texlive/texmf-dist/fonts/ty
+pe1/public/amsfonts/cm/cmr5.pfb></usr/share/texlive/texmf-dist/fonts/type1/publ
+ic/amsfonts/cm/cmr7.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfo
+nts/cm/cmsy10.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm
+/cmsy5.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmsy7.
+pfb>
+Output written on explicit.pdf (1 page, 82543 bytes).
+PDF statistics:
+ 48 PDF objects out of 1000 (max. 8388607)
+ 34 compressed objects within 1 object stream
+ 0 named destinations out of 1000 (max. 500000)
+ 1 words of extra memory for PDF output out of 10000 (max. 10000000)
+

src/time-stepping/state/sliplawstateupdater.cc → dune/tectonic/time-stepping/state/sliplawstateupdater.cc

 #include <cmath>
 
 #include "sliplawstateupdater.hh"
-#include "../../tobool.hh"
+#include "../../utils/tobool.hh"
 
 template <class ScalarVector, class Vector>
 SlipLawStateUpdater<ScalarVector, Vector>::SlipLawStateUpdater(
-    ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes,
-    double _L, double _V0)
-    : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
+        const ScalarVector& alpha_initial,
+        const BitVector& nodes,
+        const double L,
+        const double V0) :
+    nodes_(nodes),
+    L_(L),
+    V0_(V0) {
+
+    localToGlobal_.resize(nodes_.count());
+    alpha_.resize(localToGlobal_.size());
+
+    size_t localIdx = 0;
+    for (size_t i=0; i<nodes_.size(); i++) {
+        if (not toBool(nodes_[i]))
+            continue;
+
+        localToGlobal_[localIdx] = i;
+        alpha_[localIdx] = alpha_initial[i];
+        localIdx++;
+    }
+}
 
 template <class ScalarVector, class Vector>
 void SlipLawStateUpdater<ScalarVector, Vector>::nextTimeStep() {
-  alpha_o = alpha;
+    alpha_o_ = alpha_;
 }
 
 template <class ScalarVector, class Vector>
-void SlipLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
-  tau = _tau;
+void SlipLawStateUpdater<ScalarVector, Vector>::setup(double tau) {
+    tau_ = tau;
 }
 
 template <class ScalarVector, class Vector>
 void SlipLawStateUpdater<ScalarVector, Vector>::solve(
-    Vector const &velocity_field) {
-  for (size_t i = 0; i < nodes.size(); ++i) {
-    if (not toBool(nodes[i]))
-      continue;
-
-    double const V = velocity_field[i].two_norm();
-    double const mtVoL = -tau * V / L;
-    alpha[i] = (V <= 0) ? alpha_o[i] : std::expm1(mtVoL) * std::log(V / V0) +
-                                           alpha_o[i] * std::exp(mtVoL);
-  }
+    const Vector& velocity_field) {
+
+    for (size_t i=0; i<alpha_.size(); ++i) {
+        auto tangentVelocity = velocity_field[localToGlobal_[i]];
+        tangentVelocity[0] = 0.0;
+
+        double const V = tangentVelocity.two_norm();
+        double const mtVoL = -tau_ * V / L_;
+        alpha_[i] = (V <= 0) ? alpha_o_[i] : std::expm1(mtVoL) * std::log(V / V0_) +
+                                               alpha_o_[i] * std::exp(mtVoL);
+    }
 }
 
 template <class ScalarVector, class Vector>
 void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha(
-    ScalarVector &_alpha) {
-  _alpha = alpha;
+        ScalarVector& alpha) {
+
+    assert(alpha.size() == nodes_.size());
+
+    for (size_t i=0; i<localToGlobal_.size(); i++) {
+        alpha[localToGlobal_[i]] = alpha_[i];
+    }
 }
 
 template <class ScalarVector, class Vector>
-std::shared_ptr<StateUpdater<ScalarVector, Vector>>
-SlipLawStateUpdater<ScalarVector, Vector>::clone() const {
+auto SlipLawStateUpdater<ScalarVector, Vector>::clone() const
+-> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> {
   return std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(*this);
-}
+}                                                      

dune/tectonic/time-stepping/state/sliplawstateupdater.hh

+#ifndef SRC_TIME_STEPPING_STATE_SLIPLAWSTATEUPDATER_HH
+#define SRC_TIME_STEPPING_STATE_SLIPLAWSTATEUPDATER_HH
+
+#include <dune/common/bitsetvector.hh>
+
+#include "stateupdater.hh"
+
+template <class ScalarVector, class Vector>
+class SlipLawStateUpdater : public LocalStateUpdater<ScalarVector, Vector> {
+private:
+    using BitVector = Dune::BitSetVector<1>;
+
+public:
+    SlipLawStateUpdater(
+            const ScalarVector& alpha_initial,
+            const BitVector& nodes,
+            const double L,
+            const double V0);
+
+  void nextTimeStep() override;
+  void setup(double tau) override;
+  void solve(const Vector& velocity_field) override;
+  void extractAlpha(ScalarVector&) override;
+
+  auto clone() const -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> override;
+
+private:
+  std::vector<int> localToGlobal_;
+
+  ScalarVector alpha_o_;
+  ScalarVector alpha_;
+  const BitVector& nodes_;
+  const double L_;
+  const double V0_;
+  double tau_;
+};
+
+#endif

dune/tectonic/time-stepping/state/stateupdater.hh

+#ifndef SRC_TIME_STEPPING_STATE_STATEUPDATER_HH
+#define SRC_TIME_STEPPING_STATE_STATEUPDATER_HH
+
+#include <memory>
+#include <vector>
+
+// state updater for each coupling
+template <class ScalarVectorTEMPLATE, class Vector> class LocalStateUpdater {
+public:
+  using ScalarVector = ScalarVectorTEMPLATE;
+
+  void setBodyIndex(const size_t bodyIdx) {
+      bodyIdx_ = bodyIdx;
+  }
+
+  auto bodyIndex() const {
+      return bodyIdx_;
+  }
+
+  void virtual nextTimeStep() = 0;
+  void virtual setup(double _tau) = 0;
+  void virtual solve(const Vector& velocity_field) = 0;
+  void virtual extractAlpha(ScalarVector& alpha) = 0;
+
+  std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> virtual clone() const = 0;
+
+private:
+  size_t bodyIdx_;
+};
+
+
+template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
+public:
+  using ScalarVector = ScalarVectorTEMPLATE;
+  using LocalStateUpdater = LocalStateUpdater<ScalarVector, Vector>;
+
+  StateUpdater(const std::vector<size_t>& leafVertexCounts) :
+    leafVertexCounts_(leafVertexCounts) {}
+
+  void addLocalUpdater(std::shared_ptr<LocalStateUpdater> localStateUpdater) {
+    localStateUpdaters_.emplace_back(localStateUpdater);
+  }
+
+  void nextTimeStep() {
+    for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+        localStateUpdaters_[i]->nextTimeStep();
+    }
+  }
+
+  void setup(double tau) {
+    for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+        localStateUpdaters_[i]->setup(tau);
+    }
+  }
+
+  void solve(const std::vector<Vector>& velocity_field) {
+    for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+        auto& localStateUpdater = localStateUpdaters_[i];
+        localStateUpdater->solve(velocity_field[localStateUpdater->bodyIndex()]);
+    }
+  }
+
+  void extractAlpha(std::vector<ScalarVector>& alpha) {
+    alpha.resize(leafVertexCounts_.size());
+    for (size_t i=0; i<alpha.size(); i++) {
+        alpha[i].resize(leafVertexCounts_[i], 0.0);
+    }
+
+    for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+        auto& localStateUpdater = localStateUpdaters_[i];
+        localStateUpdater->extractAlpha(alpha[localStateUpdater->bodyIndex()]);
+    }
+  }
+
+  std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const {
+      auto updater = std::make_shared<StateUpdater<ScalarVector, Vector>>(leafVertexCounts_);
+
+      for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+          auto localUpdater = localStateUpdaters_[i]->clone();
+          updater->addLocalUpdater(localUpdater);
+      }
+      return updater; // std::make_shared<StateUpdater<ScalarVector, Vector>>(*this);
+  }
+
+private:
+  std::vector<size_t> leafVertexCounts_;
+  std::vector<std::shared_ptr<LocalStateUpdater>> localStateUpdaters_;
+};
+#endif

dune/tectonic/time-stepping/state_tmpl.cc

+#include "../explicitvectors.hh"
+#include "../explicitgrid.hh"
+
+#include <dune/common/promotiontraits.hh>
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+
+#include "../data-structures/friction/frictioncouplingpair.hh"
+#include "../spatial-solving/contact/dualmortarcoupling.hh"
+
+using field_type = typename Dune::PromotionTraits<typename Vector::field_type,
+                                            typename DeformedGrid::ctype>::PromotedType;
+
+using MyContactCoupling = DualMortarCoupling<field_type, DeformedGrid>;
+using MyFrictionCouplingPair = FrictionCouplingPair<DeformedGrid, LocalVector, field_type>;
+
+template std::shared_ptr<StateUpdater<ScalarVector, Vector>>
+initStateUpdater<ScalarVector, Vector, MyContactCoupling, MyFrictionCouplingPair>(
+    Config::stateModel model,
+    const std::vector<ScalarVector>& alpha_initial,
+    const std::vector<std::shared_ptr<MyContactCoupling>>& contactCouplings,
+    const std::vector<std::shared_ptr<MyFrictionCouplingPair>>& couplings);

dune/tectonic/time-stepping/step.hh

+#ifndef DUNE_TECTONIC_TIME_STEPPING_STEP_HH
+#define DUNE_TECTONIC_TIME_STEPPING_STEP_HH
+
+
+#include <future>
+#include <thread>
+#include <chrono>
+#include <functional>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "../spatial-solving/makelinearsolver.hh"
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+
+#include "stepbase.hh"
+#include "adaptivetimestepper.hh"
+
+template<class IterationStepType, class NormType, class ReductionFactorContainer>
+Dune::Solvers::Criterion reductionFactorCriterion(
+      IterationStepType& iterationStep,
+      const NormType& norm,
+      ReductionFactorContainer& reductionFactors)
+{
+  double normOfOldCorrection = 1;
+  auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+
+  return Dune::Solvers::Criterion(
+      [&, lastIterate, normOfOldCorrection] () mutable {
+        double normOfCorrection = norm.diff(*lastIterate, *iterationStep.getIterate());
+        double convRate = (normOfOldCorrection > 0) ? normOfCorrection / normOfOldCorrection : 0.0;
+
+        if (convRate>1.0)
+            std::cout << "Solver convergence rate of " << convRate << std::endl;
+
+        normOfOldCorrection = normOfCorrection;
+        *lastIterate = *iterationStep.getIterate();
+
+        reductionFactors.push_back(convRate);
+        return std::make_tuple(convRate < 0, Dune::formatString(" % '.5f", convRate));
+      },
+      " reductionFactor");
+}
+
+
+template<class IterationStepType, class Functional, class ReductionFactorContainer>
+Dune::Solvers::Criterion energyCriterion(
+      const IterationStepType& iterationStep,
+      const Functional& f,
+      ReductionFactorContainer& reductionFactors)
+{
+  double normOfOldCorrection = 1;
+  auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+
+  return Dune::Solvers::Criterion(
+      [&, lastIterate, normOfOldCorrection] () mutable {
+        double normOfCorrection = std::abs(f(*lastIterate) - f(*iterationStep.getIterate())); //norm.diff(*lastIterate, *iterationStep.getIterate());
+
+        double convRate = (normOfOldCorrection != 0.0) ? 1.0 - (normOfCorrection / normOfOldCorrection) : 0.0;
+
+        if (convRate>1.0)
+            std::cout << "Solver convergence rate of " << convRate << std::endl;
+
+        normOfOldCorrection = normOfCorrection;
+        *lastIterate = *iterationStep.getIterate();
+
+        reductionFactors.push_back(convRate);
+        return std::make_tuple(convRate < 0, Dune::formatString(" % '.5f", convRate));
+      },
+      " reductionFactor");
+}
+
+template <class ReductionFactorContainer>
+void updateReductionFactors(ReductionFactorContainer& reductionFactors) {
+    const size_t s = reductionFactors.size();
+
+    //print(reductionFactors, "reduction factors: ");
+
+    if (s>allReductionFactors.size()) {
+        allReductionFactors.resize(s);
+    }
+
+    for (size_t i=0; i<reductionFactors.size(); i++) {
+        allReductionFactors[i].push_back(reductionFactors[i]);
+    }
+
+    reductionFactors.clear();
+}
+
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+class Step : protected StepBase<Factory, ContactNetwork, Updaters, ErrorNorms> {
+public:
+    enum Mode { sameThread, newThread };
+
+private:
+    using Base = StepBase<Factory, ContactNetwork, Updaters, ErrorNorms>;
+    using NBodyAssembler = typename ContactNetwork::NBodyAssembler;
+    using UpdatersWithCount = UpdatersWithCount<Updaters>;
+
+    const Updaters& oldUpdaters_;
+    const NBodyAssembler& nBodyAssembler_;
+
+    double relativeTime_;
+    double relativeTau_;
+
+    IterationRegister& iterationRegister_;
+
+    std::packaged_task<UpdatersWithCount()> task_;
+    std::future<UpdatersWithCount> future_;
+    std::thread thread_;
+
+    Mode mode_;
+
+public:
+    Step(const Base& stepFactory, const Updaters& oldUpdaters, const NBodyAssembler& nBodyAssembler, double rTime, double rTau, IterationRegister& iterationRegister) :
+        Base(stepFactory.parset_, stepFactory.contactNetwork_, stepFactory.ignoreNodes_, stepFactory.globalFriction_,
+             stepFactory.bodywiseNonmortarBoundaries_, stepFactory.externalForces_, stepFactory.errorNorms_),
+        oldUpdaters_(oldUpdaters),
+        nBodyAssembler_(nBodyAssembler),
+        relativeTime_(rTime),
+        relativeTau_(rTau),
+        iterationRegister_(iterationRegister) {}
+
+    UpdatersWithCount doStep() {
+        // make linear solver for linear correction in TNNMGStep
+        using Vector = typename Factory::Vector;
+        using Matrix = typename Factory::Matrix;
+
+        auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(this->parset_, this->contactNetwork_);
+
+        Vector x;
+        x.resize(nBodyAssembler_.totalHasObstacle_.size());
+        x = 0;
+
+        linearSolver->getIterationStep().setProblem(x);
+      //dynamic_cast<Dune::Solvers::IterationStep<Vector>*>(linearMultigridStep.get())->setProblem(x);
+      //dynamic_cast<Dune::Solvers::IterationStep<Vector>*>(cgStep.get())->setProblem(x);
+
+      //std::vector<double> reductionFactors;
+      //linearSolver->addCriterion(reductionFactorCriterion(linearSolver->getIterationStep(), linearSolver->getErrorNorm(), reductionFactors));
+      //linearSolver->addCriterion(reductionFactorCriterion(*cgStep, norm, reductionFactors));
+
+      UpdatersWithCount newUpdatersAndCount = {oldUpdaters_.clone(), {}};
+
+      typename Base::MyCoupledTimeStepper coupledTimeStepper(
+        this->finalTime_, this->parset_, nBodyAssembler_,
+        this->ignoreNodes_, this->globalFriction_, this->bodywiseNonmortarBoundaries_,
+        newUpdatersAndCount.updaters, this->errorNorms_, this->externalForces_);
+
+      newUpdatersAndCount.count = coupledTimeStepper.step(linearSolver, relativeTime_, relativeTau_);
+      iterationRegister_.registerCount(newUpdatersAndCount.count);
+
+      //updateReductionFactors(reductionFactors);
+
+      return newUpdatersAndCount;
+    }
+
+   /* auto simple = [&]{
+        std::cout << "starting task... " << std::endl;
+        UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
+        std::this_thread::sleep_for(std::chrono::milliseconds(10000));
+        std::cout << "finishing task... " << std::endl;
+        return newUpdatersAndCount;
+    }*/
+
+    void run(Mode mode = sameThread) {
+        mode_ = mode;
+        task_ = std::packaged_task<UpdatersWithCount()>( [this]{ return doStep(); });
+        future_ = task_.get_future();
+
+        if (mode == Mode::sameThread) {
+            task_();
+        } else {
+            thread_ = std::thread(std::move(task_));
+        }
+    }
+
+    auto get() {
+        //std::cout << "Step::get() called" << std::endl;
+        if (thread_.joinable()) {
+            thread_.join();
+            //std::cout << "Thread joined " << std::endl;
+        }
+        return future_.get();
+    }
+
+    ~Step() {
+        if (thread_.joinable()) {
+            thread_.join();
+            std::cout << "Destructor:: Thread joined " << std::endl;
+        }
+    }
+};
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+Step<Factory, ContactNetwork, Updaters, ErrorNorms>
+buildStep(const StepBase<Factory, ContactNetwork, Updaters, ErrorNorms>& base,
+          const Updaters& oldUpdaters,
+          const typename ContactNetwork::NBodyAssembler& nBodyAssembler,
+          double rTime, double rTau) {
+
+    return Step(base, oldUpdaters, nBodyAssembler, rTime, rTau);
+}
+
+#endif

dune/tectonic/time-stepping/stepbase.hh

+#ifndef DUNE_TECTONIC_TIME_STEPPING_STEPBASE_HH
+#define DUNE_TECTONIC_TIME_STEPPING_STEPBASE_HH
+
+#include "coupledtimestepper.hh"
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+class StepBase {
+protected:
+  using NBodyAssembler = typename ContactNetwork::NBodyAssembler;
+  using IgnoreVector = typename Factory::BitVector;
+
+  using MyCoupledTimeStepper = CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>;
+
+  using GlobalFriction = typename MyCoupledTimeStepper::GlobalFriction;
+  using BitVector = typename MyCoupledTimeStepper::BitVector;
+  using ExternalForces = typename MyCoupledTimeStepper::ExternalForces;
+
+public:
+  StepBase(
+    Dune::ParameterTree const &parset,
+    ContactNetwork& contactNetwork,
+    const IgnoreVector& ignoreNodes,
+    GlobalFriction& globalFriction,
+    const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+    ExternalForces& externalForces,
+    const ErrorNorms& errorNorms) :
+        parset_(parset),
+        finalTime_(parset_.get<double>("problem.finalTime")),
+        contactNetwork_(contactNetwork),
+        ignoreNodes_(ignoreNodes),
+        globalFriction_(globalFriction),
+        bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
+        externalForces_(externalForces),
+        errorNorms_(errorNorms) {}
+
+  Dune::ParameterTree const &parset_;
+  double finalTime_;
+
+  ContactNetwork& contactNetwork_;
+  const IgnoreVector& ignoreNodes_;
+
+  GlobalFriction& globalFriction_;
+  const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
+
+  ExternalForces& externalForces_;
+  const ErrorNorms& errorNorms_;
+};
+
+#endif

dune/tectonic/time-stepping/uniformtimestepper.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+
+#include "uniformtimestepper.hh"
+#include "step.hh"
+
+
+/*
+ * Implementation: AdaptiveTimeStepper
+ */
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+UniformTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::UniformTimeStepper(
+        const StepBase& stepBase,
+        ContactNetwork& contactNetwork,
+        Updaters &current,
+        double relativeTime,
+        double relativeTau)
+    : relativeTime_(relativeTime),
+      relativeTau_(relativeTau),
+      stepBase_(stepBase),
+      contactNetwork_(contactNetwork),
+      current_(current) {}
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+bool UniformTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::reachedEnd() {
+  return relativeTime_ >= 1.0;
+}
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+IterationRegister UniformTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::advance() {
+  //std::cout << "AdaptiveTimeStepper::advance()" << std::endl;
+
+  using Step = Step<Factory, ContactNetwork, Updaters, ErrorNorms>;
+  iterationRegister_.reset();
+
+  UpdatersWithCount N;
+
+  auto step = Step(stepBase_, current_, contactNetwork_.nBodyAssembler(), relativeTime_, relativeTau_, iterationRegister_);
+  step.run(Step::Mode::sameThread);
+  N = step.get();
+
+  current_ = N.updaters;
+  iterationRegister_.registerFinalCount(N.count);
+  relativeTime_ += relativeTau_;
+
+  return iterationRegister_;
+}
+
+#include "uniformtimestepper_tmpl.cc"

dune/tectonic/time-stepping/uniformtimestepper.hh

+#ifndef SRC_TIME_STEPPING_UNIFORMTIMESTEPPER_HH
+#define SRC_TIME_STEPPING_UNIFORMTIMESTEPPER_HH
+
+#include "../spatial-solving/fixedpointiterator.hh"
+#include "adaptivetimestepper.hh"
+#include "stepbase.hh"
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+class UniformTimeStepper {
+
+  using UpdatersWithCount = UpdatersWithCount<Updaters>;
+  using StepBase = StepBase<Factory, ContactNetwork, Updaters, ErrorNorms>;
+
+public:
+  UniformTimeStepper(const StepBase& stepBase,
+                      ContactNetwork& contactNetwork,
+                      Updaters &current,
+                      double relativeTime,
+                      double relativeTau);
+
+  bool reachedEnd();
+
+  IterationRegister advance();
+
+  double relativeTime_;
+  const double relativeTau_;
+
+private:
+  const StepBase& stepBase_;
+  ContactNetwork& contactNetwork_;
+
+  Updaters &current_;
+
+  IterationRegister iterationRegister_;
+};
+
+#endif

dune/tectonic/time-stepping/uniformtimestepper_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include <future>
+#include <thread>
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "../spatial-solving/tnnmg/functional.hh"
+#include "../spatial-solving/solverfactory.hh"
+
+#include "../data-structures/network/contactnetwork.hh"
+#include "../data-structures/friction/globalfriction.hh"
+#include "../spatial-solving/tnnmg/zerononlinearity.hh"
+
+#include "rate/rateupdater.hh"
+#include "state/stateupdater.hh"
+#include "updaters.hh"
+
+using MyContactNetwork = ContactNetwork<Grid, Vector>;
+
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+
+using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
+
+using LinearSolver = Dune::Solvers::LoopSolver<Vector>;
+using MyNBodyAssembler = typename MyContactNetwork::NBodyAssembler;
+using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
+
+using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
+using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
+using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
+
+template class UniformTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
+
+
+using NoFriction = ZeroNonlinearity;
+using NoFrictionFunctional = Functional<Matrix&, Vector&, NoFriction&, Vector&, Vector&, double>;
+using NoFrictionSolverFactory = SolverFactory<NoFrictionFunctional, BitVector>;
+
+template class UniformTimeStepper<NoFrictionSolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
+
+/*
+template std::packaged_task<typename AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>::UpdatersWithCount()>
+AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>::step<LinearSolver>(
+        const MyUpdaters&, const MyNBodyAssembler&, std::shared_ptr<LinearSolver>&, double, double); */

dune/tectonic/utils/CMakeLists.txt

+add_custom_target(tectonic_dune_utils SOURCES
+  almostequal.hh
+  debugutils.hh 
+  diameter.hh
+  geocoordinate.hh
+  index-in-sorted-range.hh
+  maxeigenvalue.hh
+  tobool.hh
+)
+
+#install headers
+install(FILES
+  almostequal.hh
+  debugutils.hh 
+  diameter.hh
+  geocoordinate.hh
+  index-in-sorted-range.hh
+  maxeigenvalue.hh
+  tobool.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/utils/almostequal.hh

+#ifndef SRC_ALMOSTEQUAL_HH
+#define SRC_ALMOSTEQUAL_HH
+
+#include <type_traits>
+#include <limits>
+#include <math.h>
+
+template <typename ctype>
+typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type almost_equal(ctype x, ctype y, int ulp) {
+    return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp || std::abs(x-y) < std::numeric_limits<ctype>::min();
+}
+
+#endif

dune/tectonic/utils/debugutils.hh

+#ifndef DEBUG_UTILS_
+#define DEBUG_UTILS_
+
+#include <string>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+//using namespace std;
+
+    template <class Friction, typename VectorType>
+    void printRegularityTruncation(const Friction& friction, const VectorType& x, double truncationTolerance = 1e8) {
+        using BitVector = Dune::BitSetVector<MY_DIM>;
+
+        BitVector truncationFlags(x.size());
+        truncationFlags.unsetAll();
+
+        size_t count = 0;
+        size_t vsmall = 0;
+
+
+        for (size_t i = 0; i < x.size(); ++i) {
+          //std::cout << f_.phi().regularity(i, x[i]) << " xnorm: " << x[i].two_norm() << std::endl;
+          auto tangential_x = x[i];
+          tangential_x[0] = 0.0;
+
+          if (tangential_x.two_norm()<1e-14) {
+              vsmall++;
+          }
+
+          if (friction.regularity(i, x[i]) > truncationTolerance) {
+            count++;
+          }
+        }
+
+         std::cout << "V<1e-14: " << vsmall << " regularityTruncation: " << count << std::endl;
+    }
+
+    template <typename VectorType>
+    auto average(const VectorType& vec) {
+        using BlockType = typename VectorType::block_type;
+
+        BlockType res(0.0);
+        auto len = vec.size();
+
+        if (len>0) {
+            for (size_t i=0; i<vec.size(); i++) {
+                res += vec[i];
+            }
+            res *= 1.0/len;
+        }
+
+        return res;
+    }
+
+    template <typename VectorType>
+    auto minmax(const VectorType& vec) {
+        using BlockType = typename VectorType::block_type;
+
+        BlockType min(0.0), max(0.0);
+        auto len = vec.size();
+
+        if (len>0) {
+            for (size_t i=0; i<vec.size(); i++) {
+                if (vec[i] > max) {
+                    max = vec[i];
+                }
+                if (vec[i] < min) {
+                    min = vec[i];
+                }
+            }
+        }
+        return std::array{min, max};
+    }
+
+    template <int s>
+    void print(const Dune::BitSetVector<s>& x, std::string message){
+        std::cout << message << std::endl;
+        for (size_t i=0; i<x.size(); i++) {
+            std::cout << "block " << i << ": ";
+            for (size_t j=0; j<s; j++) {
+                 std::cout << x[i][j] << " ";
+            }
+            std::cout << std::endl;
+        }
+        std::cout << std::endl << std::endl;
+    }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::FieldVector<ctype, dim>& x, std::string message){
+       if (message != "")
+        std::cout << message << std::endl;
+
+       for (size_t i=0; i<x.size(); i++) {
+           std::cout << x[i] << ";" << std::endl;
+       }
+
+   }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::BlockVector<Dune::FieldVector<ctype, dim>>& x, std::string message){
+       std::cout << message  << " " << "size " << x.size() << std::endl;
+       for (size_t i=0; i<x.size(); i++) {
+           print(x[i], "");
+       }
+       std::cout << std::endl << std::endl;
+   }
+
+   template <int dim, typename ctype=double>
+   void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message, bool endOfLine = true){
+        if (message != "")
+            std::cout << message << std::endl;
+
+        for (size_t i=0; i<mat.N(); i++) {
+            for (size_t j=0; j<mat.M(); j++) {
+                if (mat.exists(i,j))
+                    std::cout << mat[i][j] << " ";
+                else
+                    std::cout << "0 ";
+            }
+
+            if (endOfLine)
+                std::cout << ";"<< std::endl;
+        }
+   }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message, const size_t line, bool endOfLine = true){
+        if (message != "")
+            std::cout << message << std::endl;
+
+        assert(line<mat.N());
+        for (size_t j=0; j<mat.M(); j++) {
+            if (mat.exists(line,j))
+                    std::cout << mat[line][j] << " ";
+                else
+                    std::cout << "0 ";
+            }
+
+            if (endOfLine)
+                std::cout << ";"<< std::endl;
+   }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::BCRSMatrix<Dune::FieldMatrix<ctype, dim, dim>>& mat, std::string message){
+        std::cout << message << " " << "size (" << mat.N() << ", " << mat.M() << ")" <<std::endl;
+
+        Dune::FieldMatrix<ctype, dim, dim> zeroEntry = 0;
+
+        for (size_t i=0; i<mat.N(); i++) {
+            for (size_t d=0; d<dim; d++) {
+                for (size_t j=0; j<mat.M(); j++) {
+                    if (mat.exists(i,j))
+                        print(mat[i][j], "", d, j==mat.M()-1);
+                    else
+                        print(zeroEntry, "", d, j==mat.M()-1);
+                }
+            }
+       }
+       std::cout << std::endl;
+   }
+
+   template <class T=Dune::FieldMatrix<double,1,1>>
+   void print(const Dune::Matrix<T>& mat, std::string message){
+       std::cout << message << std::endl;
+       for (size_t i=0; i<mat.N(); i++) {
+           for (size_t j=0; j<mat.M(); j++) {
+                   std::cout << mat[i][j] << " ";
+           }
+           std::cout << std::endl;
+       }
+       std::cout << std::endl;
+   }
+
+   template <class T>
+   void print(const std::vector<T>& x, std::string message){
+      std::cout << message << std::endl;
+      for (size_t i=0; i<x.size(); i++) {
+          std::cout << x[i] << std::endl;
+          //print(x[i], "");
+      }
+      std::cout << std::endl << std::endl;
+   }
+
+ /*  void print(const std::vector<Dune::FieldVector<double,1>>& x, std::string message){
+      std::cout << message << std::endl;
+      for (size_t i=0; i<x.size(); i++) {
+          std::cout << x[i] << " ";
+      }
+      std::cout << std::endl << std::endl;
+   }*/
+
+   template <class T>
+   void print(const std::set<T>& x, std::string message){
+      std::cout << message << std::endl;
+
+      for (const auto& c : x) {
+          std::cout << c << " ";
+      }
+      std::cout << std::endl << std::endl;
+   }
+/*
+   void print(const std::set<int>& x, std::string message){
+      std::cout << message << std::endl;
+      std::set<int>::iterator dofIt = x.begin();
+      std::set<int>::iterator dofEndIt = x.end();
+      for (; dofIt != dofEndIt; ++dofIt) {
+          std::cout << *dofIt << " ";
+      }
+      std::cout << std::endl << std::endl;
+   }
+
+   void step(const double stepNumber, std::string message=""){
+      std::cout << message << " Step " << stepNumber << "!" << std::endl;
+   }*/
+
+   template <class GridView, class VectorType>
+   void writeToVTK(const GridView& gridView, const VectorType& x, const std::string path, const std::string name) {
+       Dune::VTKWriter<GridView> vtkwriter(gridView);
+
+       std::ofstream lStream( "garbage.txt" );
+       std::streambuf* lBufferOld = std::cout.rdbuf();
+       std::cout.rdbuf(lStream.rdbuf());
+
+       vtkwriter.addVertexData(x,"data");
+       vtkwriter.pwrite(name, path, "");
+
+       std::cout.rdbuf( lBufferOld );
+   }
+
+   template <class GridView>
+   void writeToVTK(const GridView& gridView, const std::string path, const std::string name) {
+       Dune::VTKWriter<GridView> vtkwriter(gridView);
+
+       std::ofstream lStream( "garbage.txt" );
+       std::streambuf* lBufferOld = std::cout.rdbuf();
+       std::cout.rdbuf(lStream.rdbuf());
+
+       vtkwriter.pwrite(name, path, "");
+
+       std::cout.rdbuf( lBufferOld );
+   }
+
+   /*
+   template <class VectorType, class DGBasisType>
+   void writeToVTK(const DGBasisType& basis, const VectorType& x, const std::string path, const std::string name) {
+       Dune::VTKWriter<typename DGBasisType::GridView> vtkwriter(basis.getGridView());
+       VectorType toBePlotted(x);
+
+       toBePlotted.resize(basis.getGridView().indexSet().size(DGBasisType::GridView::Grid::dimension));
+
+       std::ofstream lStream( "garbage.txt" );
+       std::streambuf* lBufferOld = std::cout.rdbuf();
+       std::cout.rdbuf( lStream.rdbuf() );
+
+       vtkwriter.addVertexData(toBePlotted,"data");
+       vtkwriter.pwrite(name, path, path);
+
+       std::cout.rdbuf( lBufferOld );
+   }*/
+
+   template <class Intersection, class Basis>
+   void printIntersection(const Intersection& it, const Basis& basis, std::string message = "") {
+       if (message != "") {
+           std::cout << message << std::endl;
+       }
+
+       const auto& inside = it.inside();
+       const auto& localCoefficients = basis.getLocalFiniteElement(inside).localCoefficients();
+
+       std::cout << "dofs: ";
+       for (size_t i=0; i<localCoefficients.size(); i++) {
+           unsigned int entity = localCoefficients.localKey(i).subEntity();
+           unsigned int codim  = localCoefficients.localKey(i).codim();
+
+           if (it.containsInsideSubentity(entity, codim)) {
+               int dofIndex = basis.index(it.inside(), i);
+               std::cout << dofIndex << " ";
+               break;
+           }
+       }
+       std::cout << std::endl;
+   }
+
+   template <class BoundaryPatch, class Basis>
+   void printBoundary(const BoundaryPatch& patch, const Basis& basis, std::string message = "") {
+       if (message != "") {
+           std::cout << message << std::endl;
+       }
+
+       for (auto bIt = patch.begin(); bIt != patch.end(); ++bIt) {
+            printIntersection(*bIt, basis, "");
+       }
+
+       std::cout << std::endl;
+   }
+
+   template <class BasisType, typename ctype=double>
+   void printBasisDofLocation(const BasisType& basis) {
+     /*  typedef typename BasisType::GridView GridView;
+
+       const int dim = GridView::dimension;
+
+       std::map<int, int> indexTransformation;
+       std::map<int, Dune::FieldVector<ctype, dim>> indexCoords;
+
+
+       const GridView& gridView = basis.getGridView();
+       const int gridN = std::pow(gridView.indexSet().size(dim), 1.0/dim)-1;
+
+       for (auto& it : elements(gridView)) {
+           const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(it);
+           const auto& geometry = it.geometry();
+
+           for(int i=0; i<geometry.corners(); ++i) {
+               const auto& vertex = geometry.corner(i);
+               const auto& local = geometry.local(vertex);
+
+               int vertexIndex = vertex[0]*gridN;
+               for (int j=1; j<dim; ++j){
+                   vertexIndex += vertex[j]*gridN*std::pow(gridN+1, j);
+               }
+
+               const int localBasisSize = tFE.localBasis().size();
+               std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+               tFE.localBasis().evaluateFunction(local, localBasisRep);
+
+               for(int k=0; k<localBasisSize; ++k) {
+                   if (localBasisRep[k]==1) {
+                       int dofIndex = basis.index(it, k);
+
+                       if (indexTransformation.count(dofIndex)==0) {
+                           indexTransformation[dofIndex] = vertexIndex;
+                           indexCoords[dofIndex] = vertex;
+                       }
+
+                       break;
+                   }
+               }
+           }
+       }
+
+       std::cout << std::endl;
+       std::cout << "Coarse level: Geometric vertex indices vs. associated basis dof indices " << indexTransformation.size() << std::endl;
+       std::map<int,int>::iterator mapIt = indexTransformation.begin();
+       std::map<int,int>::iterator mapEndIt = indexTransformation.end();
+       for (; mapIt!=mapEndIt; ++mapIt) {
+           std::cout << mapIt->second << " => " << mapIt->first << " Coords: ";
+
+           const Dune::FieldVector<ctype, dim>& coords = indexCoords[mapIt->first];
+           for (size_t i=0; i<coords.size(); i++) {
+               std::cout << coords[i] << " ";
+           }
+           std::cout << std::endl;
+       }
+       std::cout << std::endl;*/
+
+       const int dim = BasisType::GridView::dimension;
+
+       std::map<int, Dune::FieldVector<ctype, dim>> indexCoords;
+
+
+       const auto& gridView = basis.getGridView();
+
+       for (auto& it : elements(gridView)) {
+            const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(it);
+            const auto& geometry = it.geometry();
+
+            for(int i=0; i<geometry.corners(); ++i) {
+                const auto& vertex = geometry.corner(i);
+                const auto& local = geometry.local(vertex);
+
+                const int localBasisSize = tFE.localBasis().size();
+                std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+                tFE.localBasis().evaluateFunction(local, localBasisRep);
+
+                for(int k=0; k<localBasisSize; ++k) {
+                    if (1.0 - localBasisRep[k]< 1e-14) {
+                        int dofIndex = basis.index(it, k);
+                        indexCoords[dofIndex] = vertex;
+                        break;
+                    }
+                }
+              }
+       }
+
+       std::cout << "Coords of basis dofs: " << indexCoords.size() << std::endl;
+       auto&& mapIt = indexCoords.begin();
+       const auto& mapEndIt = indexCoords.end();
+       for (; mapIt!=mapEndIt; ++mapIt) {
+            std::cout << mapIt->first << " Coords: ";
+
+            const auto& coords = mapIt->second;
+            for (size_t i=0; i<coords.size(); i++) {
+                std::cout << coords[i] << " ";
+            }
+            std::cout << std::endl;
+       }
+
+       std::cout << std::endl;
+   }
+
+   template <class GridView>
+   void printDofLocation(const GridView& gridView) {
+
+       std::cout << "-- GridView vertices --" << std::endl;
+       std::cout << "Index   Coords         " << std::endl;
+       std::cout << "-----------------------" << std::endl;
+
+       Dune::MultipleCodimMultipleGeomTypeMapper<
+           GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
+
+       for (auto&& it : vertices(gridView)) {
+            const auto i = vertexMapper.index(it);
+            const auto& geometry = it.geometry();
+            const auto& corner = geometry.corner(0);
+
+            std::cout << std::setw(5) << i << "   ";
+            for(size_t i=0; i<corner.size(); ++i) {
+                std::cout << std::setw(5) << std::setprecision(3) << corner[i] << " ";
+            }
+            std::cout << std::endl;
+       }
+       std::cout << std::endl;
+   }
+
+   template <class Vector, class NBodyAssembler>
+   void printMortarBasis(const NBodyAssembler& nBodyAssembler) {
+       std::cout << "-- Mortar basis in canonical coords --" << std::endl;
+       std::cout << "--------------------------------------" << std::endl;
+
+       const auto& dim = nBodyAssembler.getTransformationMatrix().N();
+
+       Vector mortarBasisVector(dim);
+       std::vector<Vector> canonicalBasisVector(nBodyAssembler.nGrids());
+       for (size_t i=0; i<dim; i++) {
+            mortarBasisVector = 0;
+            mortarBasisVector[i] = 1;
+
+            nBodyAssembler.postprocess(mortarBasisVector, canonicalBasisVector);
+            print(canonicalBasisVector, "canonicalBasisVector " + std::to_string(i));
+       }
+       std::cout << std::endl;
+   }
+#endif