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
+)
+
+#install headers
+install(FILES
+  bc.hh
+  gravity.hh 
+  midpoint.hh
+  mybody.hh
+  myglobalfrictiondata.hh
+  patchfunction.hh
+  segmented-function.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

src/one-body-problem-data/bc.hh → dune/tectonic/problem-data/bc.hh

 #ifndef SRC_ONE_BODY_PROBLEM_DATA_BC_HH
 #define SRC_ONE_BODY_PROBLEM_DATA_BC_HH
 
+#include <dune/common/function.hh>
+
 class VelocityDirichletCondition
     : public Dune::VirtualFunction<double, double> {
   void evaluate(double const &relativeTime, double &y) const {
     // Assumed to vanish at time zero
-    double const finalVelocity = -5e-5;
+      double const finalVelocity = 5e-5;
+    
+    //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 <= 0.1)
             ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
             : finalVelocity;
@@ -13,6 +23,11 @@ class VelocityDirichletCondition
 };
 
 class NeumannCondition : public Dune::VirtualFunction<double, double> {
-  void evaluate(double const &relativeTime, double &y) const { y = 0.0; }
+  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
+)
+
+#install headers
+install(FILES
+  cube.hh
+  cubefaces.hh
+  cubegridconstructor.hh
+  cuboidgeometry.hh
+  gridconstructor.hh
+  mygrids.hh
+  simplexmanager.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>;

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

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "cubefaces.hh"
+
+template <class GridView>
+template <class Vector>
+bool CubeFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
+                                    Vector const &c) {
+  return isClose2((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+}
+
+template <class GridView>
+template <class Vector>
+bool CubeFaces<GridView>::xyBoxed(Vector const &v1, Vector const &v2,
+                                Vector const &x) {
+  auto const minmax0 = std::minmax(v1[0], v2[0]);
+  auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+  if (minmax0.first - 1e-14 * lengthScale_ > x[0] or
+      x[0] > minmax0.second + 1e-14 * lengthScale_)
+    return false;
+  if (minmax1.first - 1e-14 * lengthScale_ > x[1] or
+      x[1] > minmax1.second + 1e-14 * lengthScale_)
+    return false;
+
+  return true;
+}
+
+template <class GridView>
+template <class Vector>
+bool CubeFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
+                                  Vector const &x) {
+  return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+}
+
+template <class GridView>
+CubeFaces<GridView>::CubeFaces(const GridView& gridView, const Cube& cube, double lengthScale)
+      :
+  #if MY_DIM == 3
+        lower(gridView),
+        right(gridView),
+        upper(gridView),
+        left(gridView),
+        front(gridView),
+        back(gridView),
+  #else
+        lower(gridView),
+        right(gridView),
+        upper(gridView),
+        left(gridView),
+  #endif
+        cube_(cube),
+        lengthScale_(lengthScale)
+  {
+    lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.A, cuboidGeometry.B, in.geometry().center());
+    });
+
+    right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.B, cuboidGeometry.C, in.geometry().center());
+    });
+
+    upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.D, cuboidGeometry.C, in.geometry().center());
+    });
+
+    left.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.A, cuboidGeometry.D, in.geometry().center());
+    });
+
+  #if MY_DIM == 3
+    front.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return isClose(cuboidGeometry.depth / 2.0, in.geometry().center()[2]);
+    });
+
+    back.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return isClose(-cuboidGeometry.depth / 2.0, in.geometry().center()[2]);
+    });
+  #endif
+  }
+
+#include "cubefaces_tmpl.cc"

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

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBEFACES_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBEFACES_HH
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "cube.hh"
+
+template <class GridView>
+class CubeFaces {
+private:
+    using Cube = Cube<GridView::dimensionworld>;
+
+    bool isClose(double a, double b) {
+      return std::abs(a - b) < 1e-14 * lengthScale_;
+    }
+
+    bool isClose2(double a, double b) {
+      return std::abs(a - b) <
+             1e-14 * lengthScale_ * lengthScale_;
+    }
+
+    template <class Vector>
+    bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x);
+
+    template <class Vector>
+    bool xyCollinear(Vector const &a, Vector const &b, Vector const &c);
+
+    template <class Vector>
+    bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x);
+
+public:
+    CubeFaces(const GridView& gridView, const Cube<GridView::dimensionworld>& cube, double lengthScale);
+
+    const BoundaryPatch<GridView>& lower() const {
+        return lower_;
+    }
+
+    const BoundaryPatch<GridView>& right() const {
+        return right_;
+    }
+
+    const BoundaryPatch<GridView>& upper() const {
+        return upper_;
+    }
+
+    const BoundaryPatch<GridView>& left() const {
+        return left_;
+    }
+
+#if MY_DIM == 3
+    const BoundaryPatch<GridView>& front() const {
+        return front_;
+    }
+    const BoundaryPatch<GridView>& back() const {
+        return back_;
+    }
+#endif
+
+private:
+  BoundaryPatch<GridView> lower_;
+  BoundaryPatch<GridView> right_;
+  BoundaryPatch<GridView> upper_;
+  BoundaryPatch<GridView> left_;
+
+#if MY_DIM == 3
+  BoundaryPatch<GridView> front_;
+  BoundaryPatch<GridView> back_;
+#endif
+
+  const Cube& cube_;
+  const double lengthScale_;
+};
+
+#endif
+
+

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

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+
+template struct CubeFaces<DeformedGrid::LeafGridView>;
+template struct CubeFaces<DeformedGrid::LevelGridView>;

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

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBEGRIDCONSTRUCTOR_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBEGRIDCONSTRUCTOR_HH
+
+#include "gridconstructor.hh"
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+
+template <class GridType>
+class CubeGridConstructor : public GridConstructor {    
+public:
+    using Cube = Cube<GridType::dimensionworld>;
+    
+    CubeGridConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries_);
+
+    template <class GridView>
+    void constructFaces(const GridView& gridView, CuboidGeometry const &cuboidGeometry, CubeFaces<GridView>& cubeFaces);
+
+private:
+  std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries;
+};
+
+
+#endif
+
+

src/one-body-problem-data/mygeometry.cc → dune/tectonic/problem-data/grid/cuboidgeometry.cc

@@ -10,27 +10,87 @@
 #include <cairomm/surface.h>
 #endif
 
-#include "mygeometry.hh"
-
-void MyGeometry::write() {
-  std::fstream writer("geometry", std::fstream::out);
-  writer << "A = " << A << std::endl;
-  writer << "B = " << B << std::endl;
-  writer << "C = " << C << std::endl;
-  writer << "Y = " << Y << std::endl;
-  writer << "X = " << X << std::endl;
-  writer << "Z = " << Z << std::endl;
-  writer << "U = " << U << std::endl;
-  writer << "K = " << K << std::endl;
-  writer << "M = " << M << std::endl;
-  writer << "G = " << G << std::endl;
-  writer << "H = " << H << std::endl;
-  writer << "J = " << J << std::endl;
-  writer << "I = " << I << std::endl;
-  writer << "zenith = " << zenith << std::endl;
+#include "cuboidgeometry.hh"
+
+#if MY_DIM == 3
+template <class ctype>
+CuboidGeometry<ctype>::CuboidGeometry(const GlobalCoords& origin,
+                                      const double length = 1.00, const double height = 0.27, const double depth = 0.60) :
+    length_(length*lengthScale()),
+    height_(height*lengthScale()),
+    depth_(depth*lengthScale()),
+    lowerLeft_(origin),
+    lowerRight_({origin[0]+length_, origin[1], 0}),
+    upperRight_({origin[0]+length_, origin[1]+height_, 0}),
+    upperLeft_({origin[0], origin[1]+height_, 0})
+    {}
+#else
+
+template <class ctype>
+CuboidGeometry<ctype>::CuboidGeometry(const GlobalCoords& origin,
+                                      const double length, const double height) :
+    length_(length*lengthScale()),
+    height_(height*lengthScale()),
+    lowerLeft_(origin),
+    lowerRight_({origin[0]+length_, origin[1]}),
+    upperRight_({origin[0]+length_, origin[1]+height_}),
+    upperLeft_({origin[0], origin[1]+height_})
+    {}
+#endif
+
+template <class ctype>
+void CuboidGeometry<ctype>::addWeakeningRegion(const WeakeningRegion& weakeningRegion) {
+    weakeningRegions_.emplace_back(weakeningRegion);
+}
+
+template <class ctype>
+void CuboidGeometry<ctype>::addWeakeningPatch(const Dune::ParameterTree& parset, const GlobalCoords& vertex0, const GlobalCoords& vertex1) {
+    WeakeningRegion weakPatch;
+
+#if MY_DIM == 3 // TODO: Does not work yet
+    if (vertex0 != vertex1) {
+        weakPatch.vertices.resize(4);
+        weakPatch.vertices[0] = weakPatch.vertices[2] = vertex0;
+        weakPatch.vertices[1] = weakPatch.vertices[3] = vertex1;
+
+        for (size_t k = 0; k < 2; ++k) {
+            weakPatch.vertices[k][2] = -depth_ / 2.0;
+            weakPatch.vertices[k + 2][2] = depth_ / 2.0;
+        }
+        switch (parset.get<Config::PatchType>("patchType")) {
+            case Config::Rectangular:
+                break;
+            case Config::Trapezoidal:
+                weakPatch.vertices[1][0] += 0.05 * lengthScale();
+                weakPatch.vertices[3][0] -= 0.05 * lengthScale();
+                break;
+            default:
+                assert(false);
+        }
+        addWeakeningRegion(weakPatch);
+    }
+#else
+    if (vertex0 != vertex1) {
+        weakPatch.vertices.resize(2);
+        weakPatch.vertices[0] = vertex0;
+        weakPatch.vertices[1] = vertex1;
+        addWeakeningRegion(weakPatch);
+    }
+#endif
 }
 
-void MyGeometry::render() {
+template <class ctype>
+void CuboidGeometry<ctype>::write() const {
+    std::fstream writer("geometry", std::fstream::out);
+    writer << "lowerLeft = " << lowerLeft_ << std::endl;
+    writer << "lowerRight = " << lowerRight_ << std::endl;
+    writer << "upperRight = " << upperRight_ << std::endl;
+    writer << "upperLeft = " << upperLeft_ << std::endl;
+}
+
+/*
+template <class ctype>
+void CuboidGeometry<ctype>::render() const {
 #ifdef HAVE_CAIROMM
   std::string const filename = "geometry.png";
   double const width = 600;
@@ -137,22 +197,19 @@ void MyGeometry::render() {
     cr->set_font_face(font);
     cr->set_font_size(0.03);
 
-    label(reference::A, "A");
-    label(reference::B, "B");
-    label(reference::C, "C");
-    label(reference::K, "K");
-    label(reference::M, "M");
-    label(reference::U, "U");
-    label(reference::X, "X");
-    label(reference::Y, "Y");
-    label(reference::Z, "Z");
-    label(reference::G, "G");
-    label(reference::H, "H");
-    label(reference::J, "J");
-    label(reference::I, "I");
+    label(A, "A");
+    label(B, "B");
+    label(C, "C");
+    label(D, "D");
+    label(X, "X");
+    label(Y, "Y");
+
     cr->restore();
   }
 
   surface->write_to_png(filename);
 #endif
 }
+*/
+
+#include "cuboidgeometry_tmpl.cc"

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

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+template <class ctype = double>
+class CuboidGeometry {
+public:
+    typedef Dune::FieldVector<ctype, MY_DIM> GlobalCoords;
+    using WeakeningRegion = ConvexPolyhedron<GlobalCoords>;
+
+    static constexpr double lengthScale() {
+        return 1.0;
+    } // scaling factor
+
+private:
+    const ctype length_;
+    const ctype height_;
+#if MY_DIM == 3
+    const ctype depth_;
+#endif
+
+    // corners of cube
+    const GlobalCoords lowerLeft_;
+    const GlobalCoords lowerRight_;
+    const GlobalCoords upperRight_;
+    const GlobalCoords upperLeft_;
+
+    // weakening regions
+    std::vector<WeakeningRegion> weakeningRegions_;
+
+public:
+#if MY_DIM == 3
+    CuboidGeometry(const GlobalCoords& origin,
+                   const double length = 1.00, const double height = 0.27, const double depth = 0.60);
+
+    const auto& depth() const {
+        return depth_;
+    }
+#else
+    CuboidGeometry(const GlobalCoords& origin,
+                       const double length = 1.00, const double height = 0.27);
+#endif
+
+    void addWeakeningRegion(const WeakeningRegion& weakeningRegion);
+    void addWeakeningPatch(const Dune::ParameterTree& parset, const GlobalCoords& vertex0, const GlobalCoords& vertex1);
+
+    const auto& lowerLeft() const {
+        return lowerLeft_;
+    }
+
+    const auto& lowerRight() const {
+        return lowerRight_;
+    }
+
+    const auto& upperRight() const {
+        return upperRight_;
+    }
+
+    const auto& upperLeft() const {
+        return upperLeft_;
+    }
+
+    const auto& weakeningRegions() const {
+        return weakeningRegions_;
+    }
+
+    void write() const;
+
+   // void render() const;
+};
+#endif

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

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "cuboidgeometry.hh"
+
+template class CuboidGeometry<typename Grid::ctype>;

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

+#ifndef SRC_GRIDCONSTRUCTOR_HH
+#define SRC_GRIDCONSTRUCTOR_HH
+
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+#include "../utils/diameter.hh"
+
+template <class GridType>
+class GridConstructor {
+    private:
+        double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale) {
+            return (distance / 0.0125 / lengthScale + 1.0) * smallestDiameter;
+        }
+
+    public:
+        GridConstructor() {}
+
+        std::shared_ptr<GridType> grid() {
+            return grid_;
+        }
+
+        template <class LocalVector>
+        void refine(const ConvexPolyhedron<LocalVector>& weakPatch, double smallestDiameter, double lengthScale) {
+            bool needRefine = true;
+            while (true) {
+              needRefine = false;
+              for (auto &&e : elements(grid_->leafGridView())) {
+                auto const geometry = e.geometry();
+
+                auto const weakeningRegionDistance =
+                    distance(weakPatch, geometry, 1e-6 * lengthScale);
+                auto const admissibleDiameter =
+                    computeAdmissibleDiameter(weakeningRegionDistance, smallestDiameter, lengthScale);
+
+                if (diameter(geometry) <= admissibleDiameter)
+                  continue;
+
+                needRefine = true;
+                grid_->mark(1, e);
+              }
+              if (!needRefine)
+                break;
+
+              grid_->preAdapt();
+              grid_->adapt();
+              grid_->postAdapt();
+            }
+        }
+
+        virtual void createGrid() = 0;
+
+    private:
+        Dune::GridFactory<GridType> gridFactory_;
+        std::shared_ptr<GridType> grid_;
+};
+
+
+
+
+#endif
+
+

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

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+#include "mygrids.hh"
+#include "../midpoint.hh"
+#include "../../utils/diameter.hh"
+
+#include "simplexmanager.hh"
+
+// Fix: 3D case (still Elias code)
+template <class Grid> GridsConstructor<Grid>::GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_) :
+  cuboidGeometries(cuboidGeometries_)
+{
+  size_t const gridCount = cuboidGeometries.size();
+  grids.resize(gridCount);
+  gridFactories.resize(gridCount);
+
+  for (size_t idx=0; idx<grids.size(); idx++) {
+    const auto& cuboidGeometry = *cuboidGeometries[idx];
+
+    const auto& A = cuboidGeometry.lowerLeft();
+    const auto& B = cuboidGeometry.lowerRight();
+    const auto& C = cuboidGeometry.upperRight();
+    const auto& D = cuboidGeometry.upperLeft();
+
+    unsigned int const vc = 4;
+
+#if MY_DIM == 3
+    Dune::FieldMatrix<double, 2 * vc, MY_DIM> vertices;
+#else
+    Dune::FieldMatrix<double, vc, MY_DIM> vertices;
+#endif
+    for (size_t i = 0; i < 2; ++i) {
+#if MY_DIM == 3
+      size_t numXYplanes = 2;
+#else
+      size_t numXYplanes = 1;
+#endif
+      size_t k = 0;
+      for (size_t j = 1; j <= numXYplanes; ++j) {
+        vertices[k++][i] = A[i];
+        vertices[k++][i] = B[i];
+        vertices[k++][i] = C[i];
+        vertices[k++][i] = D[i];
+        assert(k == j * vc);
+      }
+    }
+
+#if MY_DIM == 3
+    for (size_t k = 0; k < vc; ++k) {
+      vertices[k][2] = -cuboidGeometry.depth() / 2.0;
+      vertices[k + vc][2] = cuboidGeometry.depth() / 2.0;
+    }
+#endif
+
+    for (size_t i = 0; i < vertices.N(); ++i)
+      gridFactories[idx].insertVertex(vertices[i]);
+
+    Dune::GeometryType cell;
+#if MY_DIM == 3
+    cell = Dune::GeometryTypes::tetrahedron;
+#else
+    cell = Dune::GeometryTypes::triangle;
+#endif
+
+#if MY_DIM == 3
+    SimplexManager sm(vc);
+#else
+    SimplexManager sm;
+#endif
+    sm.addFromVerticesFFB(1, 2, 0);
+    sm.addFromVerticesFFB(2, 3, 0);
+    auto const &simplices = sm.getSimplices();
+
+    // sanity-check choices of simplices
+    for (size_t i = 0; i < simplices.size(); ++i) {
+      Dune::FieldMatrix<double, MY_DIM, MY_DIM> check;
+      for (size_t j = 0; j < MY_DIM; ++j)
+        check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
+      assert(check.determinant() > 0);
+      gridFactories[idx].insertElement(cell, simplices[i]);
+    }
+
+    grids[idx] = std::shared_ptr<Grid>(gridFactories[idx].createGrid());
+    grids[idx]->setRefinementType(Grid::RefinementType::COPY);
+  }
+}
+
+template <class Grid>
+std::vector<std::shared_ptr<Grid>>& GridsConstructor<Grid>::getGrids() {
+  return grids;
+}
+
+template <class Grid>
+template <class GridView>
+MyFaces<GridView> GridsConstructor<Grid>::constructFaces(
+    GridView const &gridView, const CuboidGeometry<typename Grid::ctype>& cuboidGeometry) {
+  return MyFaces<GridView>(gridView, cuboidGeometry);
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
+                                    Vector const &c) {
+  return isClose2((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyBoxed(Vector const &v1, Vector const &v2,
+                                Vector const &x) {
+  auto const minmax0 = std::minmax(v1[0], v2[0]);
+  auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+  if (minmax0.first - 1e-14 * cuboidGeometry.lengthScale() > x[0] or
+      x[0] > minmax0.second + 1e-14 * cuboidGeometry.lengthScale())
+    return false;
+  if (minmax1.first - 1e-14 * cuboidGeometry.lengthScale() > x[1] or
+      x[1] > minmax1.second + 1e-14 * cuboidGeometry.lengthScale())
+    return false;
+
+  return true;
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
+                                  Vector const &x) {
+  return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+}
+
+template <class GridView>
+MyFaces<GridView>::MyFaces(GridView const &gridView, const CuboidGeometry<typename GridView::ctype>& cuboidGeometry_)
+      :
+  #if MY_DIM == 3
+        lower(gridView),
+        right(gridView),
+        upper(gridView),
+        left(gridView),
+        front(gridView),
+        back(gridView),
+  #else
+        lower(gridView),
+        right(gridView),
+        upper(gridView),
+        left(gridView),
+  #endif
+        cuboidGeometry(cuboidGeometry_)
+  {
+    lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.lowerLeft(), cuboidGeometry.lowerRight(), in.geometry().center());
+    });
+
+    right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.lowerRight(), cuboidGeometry.upperRight(), in.geometry().center());
+    });
+
+    upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.upperLeft(), cuboidGeometry.upperRight(), in.geometry().center());
+    });
+
+    left.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.lowerLeft(), cuboidGeometry.upperLeft(), in.geometry().center());
+    });
+
+  #if MY_DIM == 3
+    front.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return isClose(cuboidGeometry.depth() / 2.0, in.geometry().center()[2]);
+    });
+
+    back.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return isClose(-cuboidGeometry.depth() / 2.0, in.geometry().center()[2]);
+    });
+  #endif
+  }
+
+double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale) {
+  return (distance / 0.0125 / lengthScale + 1.0) * smallestDiameter;
+}
+
+template <class Grid, class LocalVector>
+void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+            double smallestDiameter, double lengthScale) {
+  bool needRefine = true;
+  while (true) {
+    needRefine = false;
+    for (auto &&e : elements(grid.leafGridView())) {
+      auto const geometry = e.geometry();
+
+      auto const weakeningRegionDistance =
+          distance(weakPatch, geometry, 1e-6 * lengthScale);
+      auto const admissibleDiameter =
+          computeAdmissibleDiameter(weakeningRegionDistance, smallestDiameter, lengthScale);
+
+      if (diameter(geometry) <= admissibleDiameter)
+        continue;
+
+      needRefine = true;
+      grid.mark(1, e);
+    }
+    if (!needRefine)
+      break;
+
+    grid.preAdapt();
+    grid.adapt();
+    grid.postAdapt();
+  }
+}
+
+#include "mygrids_tmpl.cc"

src/one-body-problem-data/mygrid.hh → dune/tectonic/problem-data/grid/mygrids.hh

-#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYGRID_HH
-#define SRC_ONE_BODY_PROBLEM_DATA_MYGRID_HH
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYGRIDS_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_MYGRIDS_HH
 
 #include <dune/common/fmatrix.hh>
 #include <dune/grid/common/gridfactory.hh>
@@ -7,29 +7,32 @@
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/geometry/convexpolyhedron.hh>
 
-#include "mygeometry.hh"
+#include "cuboidgeometry.hh"
 
 template <class GridView> struct MyFaces {
   BoundaryPatch<GridView> lower;
   BoundaryPatch<GridView> right;
   BoundaryPatch<GridView> upper;
+  BoundaryPatch<GridView> left;
 
 #if MY_DIM == 3
   BoundaryPatch<GridView> front;
   BoundaryPatch<GridView> back;
 #endif
 
-  MyFaces(GridView const &gridView);
+  MyFaces(GridView const &gridView, const CuboidGeometry<typename GridView::ctype>& cuboidGeometry_);
 
 private:
+  const CuboidGeometry<typename GridView::ctype>& cuboidGeometry;
+
   bool isClose(double a, double b) {
-    return std::abs(a - b) < 1e-14 * MyGeometry::lengthScale;
-  };
+    return std::abs(a - b) < 1e-14 * cuboidGeometry.lengthScale();
+  }
 
   bool isClose2(double a, double b) {
     return std::abs(a - b) <
-           1e-14 * MyGeometry::lengthScale * MyGeometry::lengthScale;
-  };
+           1e-14 * cuboidGeometry.lengthScale() * cuboidGeometry.lengthScale();
+  }
 
   template <class Vector>
   bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x);
@@ -41,44 +44,27 @@ template <class GridView> struct MyFaces {
   bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x);
 };
 
-class SimplexManager {
-public:
-  using SimplexList = std::vector<std::vector<unsigned int>>;
-
-#if MY_DIM == 3
-  SimplexManager(unsigned int shift);
-#endif
-
-  void addFromVerticesFBB(unsigned int U, unsigned int V, unsigned int W);
-  void addFromVerticesFFB(unsigned int U, unsigned int V, unsigned int W);
-
-  SimplexList const &getSimplices();
-
-private:
-  SimplexList simplices_;
-
-#if MY_DIM == 3
-  unsigned int const shift_;
-#endif
-};
-
-template <class Grid> class GridConstructor {
+template <class Grid> class GridsConstructor {
 public:
-  GridConstructor();
+  GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_);
 
-  std::shared_ptr<Grid> getGrid();
+  std::vector<std::shared_ptr<Grid>>& getGrids();
 
   template <class GridView>
-  MyFaces<GridView> constructFaces(GridView const &gridView);
+  MyFaces<GridView> constructFaces(const GridView& gridView, const CuboidGeometry<typename Grid::ctype>& cuboidGeometry);
 
 private:
-  Dune::GridFactory<Grid> gridFactory;
+  const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries;
+  std::vector<Dune::GridFactory<Grid>> gridFactories;
+  std::vector<std::shared_ptr<Grid>> grids;
 };
 
-double computeAdmissibleDiameter(double distance, double smallestDiameter);
+double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale);
 
 template <class Grid, class LocalVector>
 void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
-            double smallestDiameter);
+            double smallestDiameter, double lengthScale);
 
 #endif
+
+

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

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "../../explicitvectors.hh"
+#include "cuboidgeometry.hh"
+
+template class GridsConstructor<Grid>;
+
+template struct MyFaces<DeformedGrid::LeafGridView>;
+template struct MyFaces<DeformedGrid::LevelGridView>;
+
+template MyFaces<DeformedGrid::LeafGridView> GridsConstructor<Grid>::constructFaces(
+    DeformedGrid::LeafGridView const &gridView, CuboidGeometry<typename Grid::ctype> const &CuboidGeometry_);
+
+template MyFaces<DeformedGrid::LevelGridView> GridsConstructor<Grid>::constructFaces(
+    DeformedGrid::LevelGridView const &gridView, CuboidGeometry<typename Grid::ctype> const &CuboidGeometry_);
+
+template void refine<Grid, LocalVector>(
+    Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+    double smallestDiameter, double lengthScale);