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

dune/tectonic/utils/maxeigenvalue.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
+#define DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
+
+/** \file
+ * \brief max eigenvalue computations for the FieldMatrix class
+ */
+
+#include <iostream>
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 1, 1>& matrix) {
+    return matrix[0][0];
+}
+
+/** \brief calculates the eigenvalues of a symmetric field matrix
+    \param[in]  matrix matrix eigenvalues are calculated for
+    \param[out] max eigenvalue
+*/
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 2, 2>& matrix) {
+    const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
+    const K p2 = p - matrix[1][1];
+    K q = p2 * p2 + matrix[1][0] * matrix[0][1];
+    if( q < 0 && q > -1e-14 ) q = 0;
+    if (q < 0) {
+        std::cout << matrix << std::endl;
+        // Complex eigenvalues are either caused by non-symmetric matrices or by round-off errors
+        DUNE_THROW(Dune::MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
+    }
+
+    // get square root
+    q = std::sqrt(q);
+
+    return p + q;
+}
+
+/** \brief Calculates the eigenvalues of a symmetric 3x3 field matrix
+    \param[in]  matrix eigenvalues are calculated for
+    \param[out] max eigenvalue
+
+    \note If the input matrix is not symmetric the behavior of this method is undefined.
+ */
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 3, 3>& matrix) {
+    Dune::FieldVector<K, 3> eigenvalues;
+    Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
+    return eigenvalues[0];
+}
+
+  /** @} end documentation */
+
+#endif

dune/tectonic/utils/reductionfactors.hh

+extern std::vector<std::vector<double>> allReductionFactors;

src/tobool.hh → dune/tectonic/utils/tobool.hh

@@ -4,7 +4,7 @@
 #include <dune/common/bitsetvector.hh>
 
 template <class Alloc>
-bool toBool(Dune::BitSetVectorConstReference<1, Alloc> x) {
+bool toBool(const Dune::BitSetVectorConstReference<1, Alloc> x) {
   return x[0];
 }
 #endif

src/CMakeLists.txt

-set(SW_SOURCE_FILES
-  assemblers.cc
-  enumparser.cc
-  hdf5/frictionalboundary-writer.cc
-  hdf5/iteration-writer.cc
-  hdf5/patchinfo-writer.cc
-  hdf5/restart-io.cc
-  hdf5/surface-writer.cc
-  hdf5/time-writer.cc
-  one-body-problem-data/mygeometry.cc
-  one-body-problem-data/mygrid.cc
-  one-body-problem.cc
-  spatial-solving/fixedpointiterator.cc
-  spatial-solving/solverfactory.cc
-  time-stepping/adaptivetimestepper.cc
-  time-stepping/coupledtimestepper.cc
-  time-stepping/rate.cc
-  time-stepping/rate/rateupdater.cc
-  time-stepping/state.cc
-  vtk.cc
-)
+add_subdirectory("foam")
+add_subdirectory("multi-body-problem")
+add_subdirectory("strikeslip")
+
 set(UGW_SOURCE_FILES
-  assemblers.cc # FIXME
-  one-body-problem-data/mygrid.cc
-  uniform-grid-writer.cc
-  vtk.cc
+  ../dune/tectonic/assemblers.cc # FIXME
+  #../dune/tectonic/io/uniform-grid-writer.cc
+  ../dune/tectonic/io/vtk.cc
+  ../dune/tectonic/problem-data/grid/mygrids.cc
 )
 
 foreach(_dim 2 3)
-  set(_sw_target one-body-problem-${_dim}D)
-  set(_ugw_target uniform-grid-writer-${_dim}D)
-
-  add_executable(${_sw_target} ${SW_SOURCE_FILES})
-  add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
-  add_dune_ug_flags(${_sw_target})
-  add_dune_ug_flags(${_ugw_target})
+  #set(_ugw_target uniform-grid-writer-${_dim}D)
 
-  add_dune_hdf5_flags(${_sw_target})
+  #add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
+  
+  #add_dune_ug_flags(${_ugw_target})
 
-  set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
-  set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  #set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
 endforeach()

src/assemblers.hh

-#ifndef SRC_ASSEMBLERS_HH
-#define SRC_ASSEMBLERS_HH
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/function.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dune/fufem/assemblers/assembler.hh>
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wsign-compare"
-#include <dune/fufem/functionspacebases/p0basis.hh>
-#pragma clang diagnostic pop
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/globalfrictiondata.hh>
-
-#include "enums.hh"
-
-template <class GridView, int dimension> class MyAssembler {
-public:
-  using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
-  using Matrix =
-      Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
-  using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
-  using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
-
-  using CellBasis = P0Basis<GridView, double>;
-  using VertexBasis = P1NodalBasis<GridView, double>;
-
-  CellBasis const cellBasis;
-  VertexBasis const vertexBasis;
-  GridView const &gridView;
-
-private:
-  using Grid = typename GridView::Grid;
-  using LocalVector = typename Vector::block_type;
-  using LocalScalarVector = typename ScalarVector::block_type;
-
-  using LocalCellBasis = typename CellBasis::LocalFiniteElement;
-  using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
-
-  Assembler<CellBasis, CellBasis> cellAssembler;
-  Assembler<VertexBasis, VertexBasis> vertexAssembler;
-
-public:
-  MyAssembler(GridView const &gridView);
-
-  void assembleFrictionalBoundaryMass(
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarMatrix &frictionalBoundaryMass) const;
-
-  void assembleMass(Dune::VirtualFunction<
-                        LocalVector, LocalScalarVector> const &densityFunction,
-                    Matrix &M) const;
-
-  void assembleElasticity(double E, double nu, Matrix &A) const;
-
-  void assembleViscosity(
-      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
-          shearViscosity,
-      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
-          bulkViscosity,
-      Matrix &C) const;
-
-  void assembleBodyForce(
-      Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-      Vector &f) const;
-
-  void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
-                       Vector &f,
-                       Dune::VirtualFunction<double, double> const &neumann,
-                       double relativeTime) const;
-
-  void assembleWeightedNormalStress(
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarVector &weightedNormalStress, double youngModulus,
-      double poissonRatio, Vector const &displacement) const;
-
-  std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
-      Config::FrictionModel frictionModel,
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      GlobalFrictionData<dimension> const &frictionInfo,
-      ScalarVector const &weightedNormalStress) const;
-
-  void assembleVonMisesStress(double youngModulus, double poissonRatio,
-                              Vector const &u, ScalarVector &stress) const;
-};
-#endif

src/explicitgrid.hh

-#ifndef SRC_EXPLICITGRID_HH
-#define SRC_EXPLICITGRID_HH
-
-#include "gridselector.hh"
-
-using GridView = Grid::LeafGridView;
-#endif