dune/tectonic/io/hdf5/restrict.hh

+#ifndef SRC_HDF5_RESTRICT_HH
+#define SRC_HDF5_RESTRICT_HH
+
+#include <cassert>
+
+#include <dune/common/bitsetvector.hh>
+
+#include "../../utils/tobool.hh"
+
+template <class Vector, class Patch>
+Vector restrictToSurface(const Vector& v1, const Patch& patch) {
+    auto const &vertices = *patch.getVertices();
+    assert(vertices.size() == v1.size());
+
+    Vector ret(vertices.count());
+    auto target = ret.begin();
+    for (size_t i = 0; i < v1.size(); ++i)
+      if (toBool(vertices[i]))
+        *(target++) = v1[i];
+    assert(target == ret.end());
+    return ret;
+}
+#endif

src/hdf5/surface-writer.cc → dune/tectonic/io/hdf5/surface-writer.hh

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
+#ifndef SRC_HDF5_SURFACE_WRITER_HH
+#define SRC_HDF5_SURFACE_WRITER_HH
 
+#include <dune/fufem/boundarypatch.hh>
 #include "restrict.hh"
-#include "surface-writer.hh"
 
-template <class ProgramState, class GridView>
-SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
-    HDF5::Grouplike &file, Vector const &vertexCoordinates, Patch const &surface)
-    : group_(file, "surface"),
+#include <dune/fufem/hdf5/file.hh>
+#include <dune/fufem/hdf5/sequenceio.hh>
+#include <dune/fufem/hdf5/singletonwriter.hh>
+
+template <class ProgramState, class GridView> class SurfaceWriter {
+  using Vector = typename ProgramState::Vector;
+  using Patch = BoundaryPatch<GridView>;
+
+public:
+  SurfaceWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
+                Patch const &surface, size_t const id) :
+      id_(id),
+      group_(file, "surface"+std::to_string(id_)),
       surface_(surface),
       surfaceDisplacementWriter_(group_, "displacement", surface.numVertices(),
                                  Vector::block_type::dimension),
       surfaceVelocityWriter_(group_, "velocity", surface.numVertices(),
                              Vector::block_type::dimension) {
+
   auto const surfaceCoordinates = restrictToSurface(vertexCoordinates, surface);
   HDF5::SingletonWriter<2> surfaceCoordinateWriter(group_, "coordinates",
                                                  surfaceCoordinates.size(),
@@ -21,15 +30,24 @@ SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
   setEntry(surfaceCoordinateWriter, surfaceCoordinates);
 }
 
-template <class ProgramState, class GridView>
-void SurfaceWriter<ProgramState, GridView>::write(
-    ProgramState const &programState) {
-  auto const surfaceDisplacements = restrictToSurface(programState.u, surface_);
-  addEntry(surfaceDisplacementWriter_, programState.timeStep,
-           surfaceDisplacements);
+  void write(ProgramState const &programState) {
 
-  auto const surfaceVelocities = restrictToSurface(programState.v, surface_);
-  addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
-}
+      auto const surfaceDisplacements = restrictToSurface(programState.u[id_], surface_);
+      addEntry(surfaceDisplacementWriter_, programState.timeStep,
+               surfaceDisplacements);
 
-#include "surface-writer_tmpl.cc"
+      auto const surfaceVelocities = restrictToSurface(programState.v[id_], surface_);
+      addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
+    }
+
+private:
+  size_t const id_;
+
+  HDF5::Group group_;
+
+  Patch const &surface_;
+
+  HDF5::SequenceIO<2> surfaceDisplacementWriter_;
+  HDF5::SequenceIO<2> surfaceVelocityWriter_;
+};
+#endif

dune/tectonic/io/hdf5/time-writer.hh

+#ifndef DUNE_TECTONIC_HDF5_TIME_WRITER_HH
+#define DUNE_TECTONIC_HDF5_TIME_WRITER_HH
+
+#include <dune/fufem/hdf5/file.hh>
+#include <dune/fufem/hdf5/sequenceio.hh>
+
+template <class ProgramState>
+class TimeWriter {
+public:
+  TimeWriter(HDF5::Grouplike &file) : file_(file),
+      relativeTimeWriter_(file_, "relativeTime"),
+      relativeTimeIncrementWriter_(file_, "relativeTimeIncrement") {}
+
+  void write(ProgramState const &programState)  {
+      addEntry(relativeTimeWriter_, programState.timeStep,
+               programState.relativeTime);
+      addEntry(relativeTimeIncrementWriter_, programState.timeStep,
+               programState.relativeTau);
+  }
+
+private:
+  HDF5::Grouplike &file_;
+
+  HDF5::SequenceIO<0> relativeTimeWriter_;
+  HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
+};
+#endif

dune/tectonic/io/io-handler.hh

+#ifndef SRC_IO_HANDLER_HH
+#define SRC_IO_HANDLER_HH
+
+#include <filesystem>
+#include <memory>
+
+#include <dune/common/parametertree.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+//#include <dune/tectonic/explicitgrid.hh>
+
+#include <dune/fufem/hdf5/file.hh>
+
+#include "../utils/geocoordinate.hh"
+
+#include "../time-stepping/adaptivetimestepper.hh"
+
+#include "hdf5/restart-io.hh"
+#include "hdf5/iteration-writer.hh"
+#include "hdf5/time-writer.hh"
+#include "hdf5-bodywriter.hh"
+
+#include "hdf5-writer.hh"
+#include "vtk.hh"
+
+template <class Assembler, class ContactNetwork, class ProgramState>
+class IOHandler {
+private:
+    using Vector = typename Assembler::Vector;
+    using LocalVector = typename Vector::block_type;
+
+    using MyVertexBasis = typename Assembler::VertexBasis;
+    using MyCellBasis = typename Assembler::CellBasis;
+    using HDF5Writer = HDF5Writer<ProgramState, MyVertexBasis, typename Assembler::GV>;
+
+    const size_t bodyCount_;
+    std::vector<size_t> nVertices_;
+
+    std::vector<Vector> vertexCoordinates_;
+    std::vector<const MyVertexBasis* > vertexBases_;
+    std::vector<const MyCellBasis* > cellBases_;
+    std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches_;
+    std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches_;
+
+    std::unique_ptr<HDF5::File> dataFile_;
+    std::unique_ptr<HDF5Writer> dataWriter_;
+
+    bool writeVTK_;
+    bool writeData_;
+
+    bool writeRestarts_;
+    size_t restartIdx_;
+    size_t restartSpacing_;
+
+    bool printProgress_;
+
+public:
+    IOHandler(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork)
+        : bodyCount_(contactNetwork.nBodies()),
+          nVertices_(bodyCount_),
+          vertexCoordinates_(bodyCount_),
+          vertexBases_(bodyCount_),
+          cellBases_(bodyCount_),
+          weakPatches_(bodyCount_),
+          writeVTK_(parset.get<bool>("vtk.write")),
+          writeData_(parset.get<bool>("data.write")),
+          writeRestarts_(parset.get<bool>("restarts.write")),
+          restartIdx_(parset.get<size_t>("restarts.first")),
+          restartSpacing_(parset.get<size_t>("restarts.spacing")),
+          printProgress_(parset.get<bool>("printProgress")) {
+
+        for (size_t i=0; i<bodyCount_; i++) {
+            nVertices_[i] = contactNetwork.body(i)->nVertices();
+        }
+
+        for (size_t i=0; i<bodyCount_; i++) {
+          const auto& body = contactNetwork.body(i);
+          vertexBases_[i] = &(body->assembler()->vertexBasis);
+          cellBases_[i] = &(body->assembler()->cellBasis);
+
+          auto& vertexCoords = vertexCoordinates_[i];
+          vertexCoords.resize(nVertices_[i]);
+
+          auto hostLeafView = body->grid()->hostGrid().leafGridView();
+          Dune::MultipleCodimMultipleGeomTypeMapper<
+              decltype(hostLeafView), Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
+          for (auto &&v : vertices(hostLeafView))
+            vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
+        }
+
+        if (writeVTK_) {
+            if (!std::filesystem::is_directory("iterates/"))
+                std::filesystem::create_directory("iterates");
+        }
+
+
+        contactNetwork.frictionPatches(frictionPatches_);
+        dataFile_ = std::make_unique<HDF5::File>("output.h5");
+        dataWriter_ = std::make_unique<HDF5Writer>(*dataFile_, vertexCoordinates_, vertexBases_,
+                              frictionPatches_);
+    }
+
+    template <class GlobalFriction>
+    void write(const ProgramState& programState, const ContactNetwork& contactNetwork, const GlobalFriction& friction, const IterationRegister& iterationCount, bool initial = false) {
+      if (writeData_) {
+        writeData(programState, contactNetwork, friction, iterationCount, initial);
+      }
+
+      if (writeVTK_) {
+        writeVTK(programState, contactNetwork);
+      }
+
+      if (writeRestarts_) {
+        writeRestarts(programState);
+      }
+
+      if (printProgress_) {
+        std::cout << "timeStep = " << std::setw(6) << programState.timeStep
+                  << ", time = " << std::setw(12) << programState.relativeTime
+                  << ", tau = " << std::setw(12) << programState.relativeTau
+                  << std::endl;
+      }
+    }
+
+    bool read(ProgramState& programState) {
+        using ADS = RestartIO<ProgramState>;
+
+        if (restartIdx_ > 0) {// automatically adjusts the time and timestep
+            auto restartFile = std::make_unique<HDF5::File>("restarts.h5", HDF5::Access::READONLY);
+            auto restartIO = std::make_unique<ADS>(*restartFile, nVertices_);
+            restartIO->read(restartIdx_, programState);
+
+            return true;
+        } else {
+            return false;
+        }
+    }
+
+private:
+    void writeVTK(const ProgramState& programState, const ContactNetwork& contactNetwork) {
+        using ScalarVector = typename Assembler::ScalarVector;
+        std::vector<ScalarVector> stress(bodyCount_);
+
+        for (size_t i=0; i<bodyCount_; i++) {
+          const auto& body = contactNetwork.body(i);
+          body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
+                                           body->data()->getPoissonRatio(),
+                                           programState.u[i], stress[i]);
+
+        }
+
+        const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases_, vertexBases_, "iterates/");
+        vtkWriter.write(programState.timeStep, programState.u, programState.v,
+                        programState.alpha, stress);
+    }
+
+    void writeRestarts(const ProgramState& programState) {
+        if (programState.timeStep % restartSpacing_ == 0) {
+          auto restartFile = std::make_unique<HDF5::File>("restarts.h5", HDF5::Access::READWRITE);
+          auto restartIO = std::make_unique<RestartIO<ProgramState>>(*restartFile, nVertices_);
+
+          restartIO->write(programState);
+          restartFile->flush();
+        }
+    }
+
+    template <class GlobalFriction>
+    void writeData(const ProgramState& programState, const ContactNetwork& contactNetwork, const GlobalFriction& friction, const IterationRegister& iterationCount, bool initial) {
+      dataWriter_->reportSolution(programState, contactNetwork, friction);
+      if (!initial) {
+        dataWriter_->reportIterations(programState, iterationCount);
+      } else {
+        dataWriter_->reportWeightedNormalStress(programState);
+      }
+
+      dataFile_->flush();
+    }
+
+};
+#endif

src/uniform-grid-writer.cc → dune/tectonic/io/uniform-grid-writer.cc

@@ -10,12 +10,15 @@
 // #include <dune/common/parametertree.hh>
 // #include <dune/common/parametertreeparser.hh>
 
-#include "assemblers.hh"
-#include "diameter.hh"
-#include "gridselector.hh"
-#include "one-body-problem-data/mygrid.hh"
+#include "../assemblers.hh"
+#include "../gridselector.hh"
+
 #include "vtk.hh"
 
+#include "../problem-data/grid/gridconstructor.hh"
+
+#include "../utils/diameter.hh"
+
 size_t const dims = MY_DIM;
 size_t const refinements = 5; // FIXME?
 

dune/tectonic/io/vtk-bodywriter.hh

+#ifndef SRC_VTK_BODYWRITER_HH
+#define SRC_VTK_BODYWRITER_HH
+
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <dune/fufem/functions/vtkbasisgridfunction.hh>
+#include "../utils/debugutils.hh"
+
+template <class VertexBasis, class CellBasis> class VTKBodyWriter {
+private:
+  CellBasis const &cellBasis_;
+  VertexBasis const &vertexBasis_;
+  std::string const prefix_;
+
+public:
+  VTKBodyWriter(CellBasis const &cellBasis, VertexBasis const &vertexBasis, std::string prefix) :
+    cellBasis_(cellBasis), vertexBasis_(vertexBasis), prefix_(prefix)
+  {}
+
+  template <class Vector, class ScalarVector>
+  void write(size_t record, Vector const &u, Vector const &v,
+             ScalarVector const &alpha, ScalarVector const &stress) const {
+      Dune::VTKWriter<typename VertexBasis::GridView> writer(
+          vertexBasis_.getGridView());
+
+      auto const displacementPointer =
+          std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
+              vertexBasis_, u, "displacement");
+      writer.addVertexData(displacementPointer);
+
+      auto const velocityPointer =
+          std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
+              vertexBasis_, v, "velocity");
+      writer.addVertexData(velocityPointer);
+
+      auto const AlphaPointer =
+          std::make_shared<VTKBasisGridFunction<VertexBasis, ScalarVector> const>(
+              vertexBasis_, alpha, "Alpha");
+      writer.addVertexData(AlphaPointer);
+
+      auto const stressPointer =
+          std::make_shared<VTKBasisGridFunction<CellBasis, ScalarVector> const>(
+              cellBasis_, stress, "stress");
+      writer.addCellData(stressPointer);
+
+      std::string const filename = prefix_ + "_" + std::to_string(record);
+      writer.write(filename.c_str(), Dune::VTK::appendedraw);
+  }
+
+  void writeGrid() const {
+      Dune::VTKWriter<typename VertexBasis::GridView> writer(
+          vertexBasis_.getGridView());
+
+      std::string const filename = prefix_ + "_grid";
+      writer.write(filename.c_str(), Dune::VTK::appendedraw);
+  }
+};
+
+#endif

dune/tectonic/io/vtk.hh

+#ifndef SRC_VTK_HH
+#define SRC_VTK_HH
+
+#include <string>
+
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+#include <dune/fufem/functions/vtkbasisgridfunction.hh>
+
+#include "../utils/debugutils.hh"
+#include "vtk-bodywriter.hh"
+
+template <class VertexBasis, class CellBasis> class MyVTKWriter {
+private:
+  std::vector<VTKBodyWriter<VertexBasis, CellBasis>* > bodyVTKWriters_;
+  std::string const prefix_;
+
+public:
+  MyVTKWriter(const std::vector<const CellBasis* >& cellBases, const std::vector<const VertexBasis* >& vertexBases,
+              std::string prefix): bodyVTKWriters_(vertexBases.size()), prefix_(prefix) {
+
+      for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+        bodyVTKWriters_[i] = new VTKBodyWriter<VertexBasis, CellBasis>(*cellBases[i], *vertexBases[i], prefix_+std::to_string(i));
+      }
+  }
+
+  ~MyVTKWriter() {
+      for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+        delete bodyVTKWriters_[i];
+      }
+  }
+
+  template <class Vector, class ScalarVector>
+  void write(size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v,
+             const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const {
+
+      for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+        bodyVTKWriters_[i]->write(record, u[i], v[i], alpha[i], stress[i]);
+      }
+    }
+
+  void writeGrids() const {
+      for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+        bodyVTKWriters_[i]->writeGrid();
+      }
+  }
+};
+
+#endif

dune/tectonic/minimisation.hh

-#ifndef DUNE_TECTONIC_MINIMISATION_HH
-#define DUNE_TECTONIC_MINIMISATION_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/stdstreams.hh>
-
-#include <dune/fufem/arithmetic.hh>
-#include <dune/fufem/interval.hh>
-#include <dune/tnnmg/problem-classes/bisection.hh>
-
-#include <dune/tectonic/mydirectionalconvexfunction.hh>
-
-// Warning: this exploits the property v*x = 0
-template <class Functional>
-double lineSearch(Functional const &J,
-                  typename Functional::LocalVector const &x,
-                  typename Functional::LocalVector const &v,
-                  Bisection const &bisection) {
-  MyDirectionalConvexFunction<typename Functional::Nonlinearity> const JRest(
-      J.alpha * v.two_norm2(), J.b * v, J.phi, x, v);
-  int count;
-  return bisection.minimize(JRest, 0.0, 0.0, count);
-}
-
-/** Minimise a quadratic problem, for which both the quadratic and the
-    nonlinear term have gradients which point in the direction of
-    their arguments */
-template <class Functional>
-void minimise(Functional const &J, typename Functional::LocalVector &x,
-              Bisection const &bisection) {
-  auto v = J.b;
-  double const vnorm = v.two_norm();
-  if (vnorm <= 0.0)
-    return;
-  v /= vnorm;
-
-  double const alpha = lineSearch(J, x, v, bisection);
-  Arithmetic::addProduct(x, alpha, v);
-}
-#endif

dune/tectonic/multi-threading/CMakeLists.txt

+
+add_custom_target(tectonic_dune_multi-threading SOURCES
+  stoppable.hh
+  task.hh 
+)
+
+#install headers
+install(FILES
+  stoppable.hh
+  task.hh 
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/multi-threading/task.hh

+#include <future>
+#include <utility>
+
+/*
+ * Class that encapsulates promise and future object,
+ * provides API to set exit signal for thread
+ */
+class Task {
+protected:
+    std::promise<void> exitSignal;
+    std::future<void> futureObj;
+
+public:
+    Task() :
+        futureObj(exitSignal.get_future())
+    {}
+
+    Task(Task&& obj) :
+        exitSignal(std::move(obj.exitSignal)), futureObj(std::move(obj.futureObj))
+    {}
+
+    Task& operator=(Task&& obj) {
+        exitSignal = std::move(obj.exitSignal);
+        futureObj = std::move(obj.futureObj);
+        return *this;
+    }
+
+    virtual void run_task() = 0;
+
+    // thread function to be executed by thread
+    void operator()() {
+        run_task();
+    }
+
+    bool stopRequested() {
+        // checks if value in future object is available
+        if (futureObj.wait_for(std::chrono::milliseconds(0)) == std::future_status::timeout)
+            return false;
+        return true;
+    }
+
+    // stop thread by setting value in promise object
+    void stop() {
+        exitSignal.set_value();
+    }
+};
+
+
+

dune/tectonic/myblockproblem.hh

-#ifndef DUNE_TECTONIC_MYBLOCKPROBLEM_HH
-#define DUNE_TECTONIC_MYBLOCKPROBLEM_HH
-
-// Based on dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/fufem/arithmetic.hh>
-#include <dune/solvers/common/interval.hh>
-#include <dune/solvers/computeenergy.hh>
-#include <dune/tnnmg/problem-classes/bisection.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineargsproblem.hh>
-
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/minimisation.hh>
-#include <dune/tectonic/mydirectionalconvexfunction.hh>
-#include <dune/tectonic/quadraticenergy.hh>
-
-/** \brief Base class for problems where each block can be solved with a
- * modified gradient method */
-template <class ConvexProblem>
-class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
-private:
-  typedef BlockNonlinearGSProblem<ConvexProblem> Base;
-
-public:
-  using typename Base::ConvexProblemType;
-  using typename Base::LocalMatrixType;
-  using typename Base::LocalVectorType;
-  using typename Base::MatrixType;
-  using typename Base::VectorType;
-
-  size_t static const block_size = ConvexProblem::block_size;
-  size_t static const coarse_block_size = block_size;
-
-  /** \brief Solves one local system */
-  class IterateObject;
-
-  struct Linearization {
-    size_t static const block_size = coarse_block_size;
-
-    using LocalMatrix = typename MyBlockProblem<ConvexProblem>::LocalMatrixType;
-    using MatrixType = Dune::BCRSMatrix<typename Linearization::LocalMatrix>;
-    using VectorType =
-        Dune::BlockVector<Dune::FieldVector<double, Linearization::block_size>>;
-    using BitVectorType = Dune::BitSetVector<Linearization::block_size>;
-
-    typename Linearization::MatrixType A;
-    typename Linearization::VectorType b;
-    typename Linearization::BitVectorType ignore;
-
-    Dune::BitSetVector<Linearization::block_size> truncation;
-  };
-
-  MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
-      : Base(parset, problem),
-        maxEigenvalues_(problem.f.size()),
-        localBisection(0.0, 1.0, 0.0, true, 0.0) {
-    for (size_t i = 0; i < problem.f.size(); ++i) {
-      LocalVectorType eigenvalues;
-      Dune::FMatrixHelp::eigenValues(problem.A[i][i], eigenvalues);
-      maxEigenvalues_[i] =
-          *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
-    }
-  }
-
-  std::string getOutput(bool header = false) const {
-    if (header) {
-      outStream.str("");
-      for (size_t j = 0; j < block_size; ++j)
-        outStream << "  trunc" << std::setw(2) << j;
-    }
-    std::string s = outStream.str();
-    outStream.str("");
-    return s;
-  }
-
-  double computeEnergy(const VectorType &v) const {
-    return 0.0; // FIXME
-    // return ::computeEnergy(problem_.A, v, problem_.f) + problem_.phi(v);
-  }
-
-  void projectCoarseCorrection(VectorType const &u,
-                               typename Linearization::VectorType const &v,
-                               VectorType &projected_v,
-                               Linearization const &linearization) const {
-    projected_v = v;
-    for (size_t i = 0; i < v.size(); ++i)
-      for (size_t j = 0; j < block_size; ++j)
-        if (linearization.truncation[i][j])
-          projected_v[i][j] = 0;
-  }
-
-  double computeDampingParameter(VectorType const &u,
-                                 VectorType const &projected_v) const {
-    VectorType v = projected_v;
-
-    double const vnorm = v.two_norm();
-    if (vnorm <= 0)
-      return 1.0;
-
-    v /= vnorm; // Rescale for numerical stability
-
-    auto const psi = restrict(problem_.A, problem_.f, u, v, problem_.phi);
-    Dune::Solvers::Interval<double> D;
-    psi.subDiff(0, D);
-    if (D[1] > 0) // NOTE: Numerical instability can actually get us here
-      return 0;
-
-    int bisectionsteps = 0;
-    Bisection const globalBisection; // NOTE: defaults
-    return globalBisection.minimize(psi, vnorm, 0.0, bisectionsteps) / vnorm;
-  }
-
-  void assembleTruncate(VectorType const &u, Linearization &linearization,
-                        Dune::BitSetVector<block_size> const &ignore) const {
-    // we can just copy the ignore information
-    linearization.ignore = ignore;
-
-    // determine truncation pattern
-    linearization.truncation.resize(u.size());
-    linearization.truncation.unsetAll();
-    for (size_t i = 0; i < u.size(); ++i) {
-      if (problem_.phi.regularity(i, u[i]) > 1e8) { // TODO: Make customisable
-        linearization.truncation[i] = true;
-        continue;
-      }
-
-      for (size_t j = 0; j < block_size; ++j)
-        if (linearization.ignore[i][j])
-          linearization.truncation[i][j] = true;
-    }
-
-    // construct sparsity pattern for linearization
-    Dune::MatrixIndexSet indices(problem_.A.N(), problem_.A.M());
-    indices.import(problem_.A);
-    problem_.phi.addHessianIndices(indices);
-
-    // construct matrix from pattern and initialize it
-    indices.exportIdx(linearization.A);
-    linearization.A = 0.0;
-
-    // compute quadratic part of hessian (linearization.A += problem_.A)
-    for (size_t i = 0; i < problem_.A.N(); ++i) {
-      auto const end = std::end(problem_.A[i]);
-      for (auto it = std::begin(problem_.A[i]); it != end; ++it)
-        linearization.A[i][it.index()] += *it;
-    }
-
-    // compute nonlinearity part of hessian
-    problem_.phi.addHessian(u, linearization.A);
-
-    // compute quadratic part of gradient
-    linearization.b.resize(u.size());
-    problem_.A.mv(u, linearization.b);
-    linearization.b -= problem_.f;
-
-    // compute nonlinearity part of gradient
-    problem_.phi.addGradient(u, linearization.b);
-
-    // -grad is needed for Newton step
-    linearization.b *= -1.0;
-
-    // apply truncation to stiffness matrix and rhs
-    for (size_t row = 0; row < linearization.A.N(); ++row) {
-      auto const col_end = std::end(linearization.A[row]);
-      for (auto col_it = std::begin(linearization.A[row]); col_it != col_end;
-           ++col_it) {
-        size_t const col = col_it.index();
-        for (size_t i = 0; i < col_it->N(); ++i) {
-          auto const blockEnd = std::end((*col_it)[i]);
-          for (auto blockIt = std::begin((*col_it)[i]); blockIt != blockEnd;
-               ++blockIt)
-            if (linearization.truncation[row][i] or
-                linearization.truncation[col][blockIt.index()])
-              *blockIt = 0.0;
-        }
-      }
-      for (size_t j = 0; j < block_size; ++j)
-        if (linearization.truncation[row][j])
-          linearization.b[row][j] = 0.0;
-    }
-    for (size_t j = 0; j < block_size; ++j)
-      outStream << std::setw(9) << linearization.truncation.countmasked(j);
-  }
-
-  /** \brief Constructs and returns an iterate object */
-  IterateObject getIterateObject() {
-    return IterateObject(localBisection, problem_, maxEigenvalues_);
-  }
-
-private:
-  std::vector<double> maxEigenvalues_;
-
-  // problem data
-  using Base::problem_;
-
-  Bisection const localBisection;
-
-  mutable std::ostringstream outStream;
-};
-
-/** \brief Solves one local system using a scalar Gauss-Seidel method */
-template <class ConvexProblem>
-class MyBlockProblem<ConvexProblem>::IterateObject {
-  friend class MyBlockProblem;
-
-protected:
-  /** \brief Constructor, protected so only friends can instantiate it
-   * \param bisection The class used to do a scalar bisection
-   * \param problem The problem including quadratic part and nonlinear part
-   */
-  IterateObject(Bisection const &bisection, ConvexProblem const &problem,
-                std::vector<double> const &maxEigenvalues)
-      : problem(problem),
-        maxEigenvalues_(maxEigenvalues),
-        bisection_(bisection) {}
-
-public:
-  /** \brief Set the current iterate */
-  void setIterate(VectorType &u) {
-    this->u = u;
-    return;
-  }
-
-  /** \brief Update the i-th block of the current iterate */
-  void updateIterate(LocalVectorType const &ui, size_t i) {
-    u[i] = ui;
-    return;
-  }
-
-  /** \brief Minimise a local problem
-   * \param[out] ui The solution
-   * \param m Block number
-   * \param ignore Set of degrees of freedom to leave untouched
-   */
-  void solveLocalProblem(
-      LocalVectorType &ui, size_t m,
-      typename Dune::BitSetVector<block_size>::const_reference ignore) {
-    {
-      LocalVectorType localb = problem.f[m];
-      auto const end = std::end(problem.A[m]);
-      for (auto it = std::begin(problem.A[m]); it != end; ++it) {
-        size_t const j = it.index();
-        Arithmetic::subtractProduct(localb, *it, u[j]); // also the diagonal!
-      }
-      Arithmetic::addProduct(localb, maxEigenvalues_[m], u[m]);
-
-      // We minimise over an affine subspace
-      for (size_t j = 0; j < block_size; ++j)
-        if (ignore[j])
-          localb[j] = 0;
-        else
-          ui[j] = 0;
-
-      QuadraticEnergy<
-          typename ConvexProblem::NonlinearityType::LocalNonlinearity>
-          localJ(maxEigenvalues_[m], localb, problem.phi.restriction(m));
-      minimise(localJ, ui, bisection_);
-    }
-  }
-
-private:
-  ConvexProblem const &problem;
-  std::vector<double> maxEigenvalues_;
-  Bisection const bisection_;
-  // state data for smoothing procedure used by:
-  // setIterate, updateIterate, solveLocalProblem
-  VectorType u;
-};
-#endif

dune/tectonic/mydirectionalconvexfunction.hh

-#ifndef DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
-#define DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
-
-// Copied from dune/tnnmg/problem-classes/directionalconvexfunction.hh
-// Allows phi to be const
-
-#include <dune/fufem/arithmetic.hh>
-#include <dune/solvers/common/interval.hh>
-
-template <class Nonlinearity> class MyDirectionalConvexFunction {
-public:
-  using Vector = typename Nonlinearity::VectorType;
-  using Matrix = typename Nonlinearity::MatrixType;
-
-  MyDirectionalConvexFunction(double A, double b, Nonlinearity const &phi,
-                              Vector const &u, Vector const &v)
-      : A(A), b(b), phi(phi), u(u), v(v) {
-    phi.directionalDomain(u, v, dom);
-  }
-
-  double quadraticPart() const { return A; }
-
-  double linearPart() const { return b; }
-
-  void subDiff(double x, Dune::Solvers::Interval<double> &D) const {
-    Vector uxv = u;
-    Arithmetic::addProduct(uxv, x, v);
-    phi.directionalSubDiff(uxv, v, D);
-    auto const Axmb = A * x - b;
-    D[0] += Axmb;
-    D[1] += Axmb;
-  }
-
-  void domain(Dune::Solvers::Interval<double> &domain) const {
-    domain[0] = this->dom[0];
-    domain[1] = this->dom[1];
-  }
-
-  double A;
-  double b;
-
-private:
-  Nonlinearity const &phi;
-  Vector const &u;
-  Vector const &v;
-
-  Dune::Solvers::Interval<double> dom;
-};
-
-/*
-  1/2 <A(u + hv),u + hv> - <b, u + hv>
-  = 1/2 <Av,v> h^2 - <b - Au, v> h + const.
-
-  localA = <Av,v>
-  localb = <b - Au, v>
-*/
-
-template <class Matrix, class Vector, class Nonlinearity>
-MyDirectionalConvexFunction<Nonlinearity> restrict(Matrix const &A,
-                                                   Vector const &b,
-                                                   Vector const &u,
-                                                   Vector const &v,
-                                                   Nonlinearity const &phi) {
-  return MyDirectionalConvexFunction<Nonlinearity>(
-      Arithmetic::Axy(A, v, v), Arithmetic::bmAxy(A, b, u, v), phi, u, v);
-}
-
-#endif

dune/tectonic/problem-data/CMakeLists.txt

+add_subdirectory("grid")
+
+add_custom_target(tectonic_dune_problem-data SOURCES
+  bc.hh
+  gravity.hh 
+  midpoint.hh
+  mybody.hh
+  myglobalfrictiondata.hh
+  patchfunction.hh
+  segmented-function.hh
+  strikeslipbody.hh
+)
+
+#install headers
+install(FILES
+  bc.hh
+  gravity.hh 
+  midpoint.hh
+  mybody.hh
+  myglobalfrictiondata.hh
+  patchfunction.hh
+  segmented-function.hh
+  strikeslipbody.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/problem-data/bc.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+
+#include <cassert>
+
+#include <dune/common/function.hh>
+
+class VelocityDirichletCondition
+    : public Dune::VirtualFunction<double, double> {
+
+public:
+  VelocityDirichletCondition(const double finalVelocity = 5e-5, const double threshold = 0.1) :
+      finalVelocity_(finalVelocity),
+      threshold_(threshold)
+  {}
+
+  void evaluate(double const &relativeTime, double &y) const {
+    // Assumed to vanish at time zero
+    
+    //std::cout << "VelocityDirichletCondition::evaluate()" << std::endl;
+    
+    /*if (relativeTime <= 0.1)
+        std::cout << "- loading phase" << std::endl;
+    else
+        std::cout << "- final velocity reached" << std::endl;*/
+    
+    y = (relativeTime <= threshold_)
+            ? finalVelocity_ * (1.0 - std::cos(relativeTime * M_PI / threshold_)) / 2.0
+            : finalVelocity_;
+
+    //y = relativeTime * finalVelocity;
+
+  }
+
+private:
+  const double finalVelocity_;
+  const double threshold_;
+};
+
+class VelocityStepDirichletCondition
+    : public Dune::VirtualFunction<double, double> {
+
+public:
+  VelocityStepDirichletCondition(const double initialVelocity = 5e-5, const double finalVelocity = 1e-4, const double loadingThreshold = 0.1, const double stepTime = 0.5) :
+      initialVelocity_(initialVelocity),
+      finalVelocity_(finalVelocity),
+      loadingThreshold_(loadingThreshold),
+      stepTime_(stepTime)
+  {
+      assert(loadingThreshold_ < stepTime_);
+  }
+
+  void evaluate(double const &relativeTime, double &y) const {
+    // Assumed to vanish at time zero
+
+    y = finalVelocity_;
+
+    if (relativeTime <= loadingThreshold_) {
+        y = initialVelocity_ * (1.0 - std::cos(relativeTime * M_PI / loadingThreshold_)) / 2.0;
+    } else if (relativeTime < stepTime_) {
+        y = initialVelocity_;
+    }
+  }
+
+private:
+  const double initialVelocity_;
+  const double finalVelocity_;
+  const double loadingThreshold_;
+  const double stepTime_;
+};
+
+class NeumannCondition : public Dune::VirtualFunction<double, double> {
+  public:
+    NeumannCondition(double c = 0.0) : c_(c) {}
+    void evaluate(double const &relativeTime, double &y) const { y = -c_; }
+
+  private:
+    double c_;
+};
+#endif

dune/tectonic/problem-data/grid/CMakeLists.txt

+add_custom_target(tectonic_dune_problem-data_grid SOURCES
+  cube.hh
+  cube.cc 
+  cubefaces.hh
+  cubefaces.cc
+  cubegridconstructor.hh
+  cuboidgeometry.hh
+  cuboidgeometry.cc
+  gridconstructor.hh
+  mygrids.hh
+  mygrids.cc
+  simplexmanager.hh
+  simplexmanager.cc
+  trianglegeometry.hh
+  trianglegeometry.cc
+  trianglegrids.hh
+  trianglegrids.cc
+)
+
+#install headers
+install(FILES
+  cube.hh
+  cubefaces.hh
+  cubegridconstructor.hh
+  cuboidgeometry.hh
+  gridconstructor.hh
+  mygrids.hh
+  simplexmanager.hh
+  trianglegeometry.hh
+  trianglegrids.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/problem-data/grid/cube.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "cube.hh"
+
+template <int dim>
+Cube<dim>::Cube(bool invariant = false) : invariant_(invariant) {
+    if (invariant_) {
+        nChildren_ = 1;
+    } else {
+        nChildren_ = std::pow(2, dim);
+    }
+    children_.resize(nChildren_, nullptr);
+}
+
+template <int dim>
+Cube<dim>::Cube(const Vector& A, const Vector& B, bool invariant = false) : invariant_(invariant) {
+    setCorners(A, B);
+
+    if (invariant_) {
+        nChildren_ = 1;
+    } else {
+        nChildren_ = std::pow(2, dim);
+    }
+    children_.resize(nChildren_, nullptr);
+}
+
+// generate a combination of unit basis vectors from binary number
+template <int dim>
+void Cube<dim>::parseBinary(int binary, Vector& shift) const {
+
+    shift = 0;
+    for (int d=0; d<dim; d++) {
+        // yields true, if d-th digit of binary is a 1
+        if (binary & std::pow(2, d)) {
+            // use d-th basis vector
+            shift[d] = 1;
+        }
+    }
+}
+
+// constructs child cubes and sets children_
+template <int dim>
+void Cube<dim>::split() {
+    if (invariant_) {
+        children_[0] = this;
+    } else {
+        const int nNewCubes = std::pow(2, dim);
+        newCubes.resize(nNewCubes);
+
+        Vector midPoint = A_;
+        midPoint += 1/2*(B_ - A_);
+
+        const double h = std::pow(2.0, std::round(std::log(midPoint[0]-A_[0])/std::log(2)));
+
+        newCubes[0] = new Cube(A_, midPoint, false);
+        newCubes[nNewCubes-1] = new Cube(midPoint, B_, true);
+
+        for (int i=1; i<nNewCubes-1; i++) {
+            Vector shift;
+            parseBinary(i, shift);
+            shift *= h;
+
+            newCubes[i] = new Cube(A_+shift, midPoint+shift);
+        }
+    } else {
+        children_.resize(1);
+    }
+}
+
+#include "cube_tmpl.cc"

dune/tectonic/problem-data/grid/cube.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBE_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBE_HH
+
+
+// works in any space dimension
+template <int dim>
+class Cube {
+private:
+    using Vector = Dune::FieldVector<double, dim>;
+
+    // generate a combination of unit basis vectors from binary number
+    void parseBinary(int binary, Vector& shift) const;
+
+public:
+    Cube(bool invariant = false);
+    Cube(const Vector& A, const Vector& B, bool invariant = false);
+
+    void setCorners(const Vector& A, const Vector& B) {
+        A_ = A;
+        B_ = B;
+    }
+
+    void setParent(Cube* parent) {
+        parent_ = parent;
+    }
+
+    const std::vector<std::shared_ptr<Cube>>& children() const {
+        return children_;
+    }
+
+    // constructs child cubes and sets children_
+    void split();
+
+private:
+    using Vector = Dune::FieldVector<double, dim>;
+
+    Vector A_; // two corners that are diagonal define dim-cube in any space dimension
+    Vector B_;
+
+    const bool invariant_; // flag to set if Cube can be split()
+
+    std::shared_ptr<Cube<dim>> parent_ = nullptr;
+    std::vector<std::shared_ptr<Cube<dim>>> children_;
+    int nChildren_;
+};
+
+#endif
+
+

dune/tectonic/problem-data/grid/cube_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+template class Cube<MY_DIM>;