.gitignore

-*.o
-*.pyc
-.clang-format
-.deps
-.libs
-/aclocal.m4
-/am
-/autom4te.cache
-/compile
-/config.guess
-/config.h
-/config.h.in
-/config.log
-/config.lt
-/config.status
-/config.sub
-/configure
-/depcomp
-/dependencies.m4
-/dune-tectonic.pc
-/install-sh
-/libtool
-/ltmain.sh
-/missing
-/stamp-h1
-/test-driver
-Makefile
-Makefile.in
-src/sand-wedge-?D
-src/sliding-block-?D
+build-cmake

CMakeLists.txt

+cmake_minimum_required(VERSION 2.8.6)
+project(dune-tectonic CXX)
+
+if(NOT (dune-common_DIR OR dune-common_ROOT OR
+      "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
+    string(REPLACE  ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR
+      ${PROJECT_BINARY_DIR})
+endif()
+
+#find dune-common and set the module path
+find_package(dune-common REQUIRED)
+list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules"
+  ${dune-common_MODULE_PATH})
+
+#include the dune macros
+include(DuneMacros)
+
+# start a dune project with information from dune.module
+dune_project()
+
+find_package(HDF5 COMPONENTS C REQUIRED)
+
+add_subdirectory("src")
+add_subdirectory("dune")
+add_subdirectory("doc")
+add_subdirectory("cmake/modules")
+
+# finalize the dune project, e.g. generating config.h etc.
+finalize_dune_project(GENERATE_CONFIG_H_CMAKE)

Makefile.am

-EXTRA_DIST = dune.module
-
-SUBDIRS = src m4 dune
-
-include $(top_srcdir)/am/top-rules
-include $(top_srcdir)/am/global-rules
-

cmake/modules/CMakeLists.txt

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

cmake/modules/DuneTectonicMacros.cmake

+# File for module specific CMake tests.

config.h.cmake

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

configure.ac

-# -*- Autoconf -*-
-# Process this file with autoconf to produce a configure script.
-AC_PREREQ(2.50)
-DUNE_AC_INIT # gets module version from dune.module file
-AM_INIT_AUTOMAKE
-AM_SILENT_RULES
-AC_CONFIG_SRCDIR([dune-tectonic.pc.in])
-AC_CONFIG_HEADERS([config.h])
-
-DUNE_CHECK_ALL
-
-AC_ARG_WITH([cairomm], AS_HELP_STRING([--with-cairomm], [Use cairomm]))
-
-AS_IF([test "x$with_cairomm" != "xno"], [
-   PKG_CHECK_MODULES([CAIROMM], [cairomm-1.0],
-                     [AM_CONDITIONAL(CAIROMM, true)
-                      AC_DEFINE([HAVE_CAIROMM], [1], [If CairoMM is available])
-   ])
-])
-
-
-AC_CONFIG_FILES([
-  Makefile
-  src/Makefile
-  dune/Makefile
-  dune/tectonic/Makefile
-  m4/Makefile
-  dune-tectonic.pc
-])
-AC_OUTPUT
-# finally print the summary information
-DUNE_SUMMARY_ALL

doc/CMakeLists.txt

+add_subdirectory("doxygen")

doc/doxygen/CMakeLists.txt

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

doc/doxygen/Doxylocal

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

dune-tectonic.pc.in

@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@
 Version: @VERSION@
 Description: dune-tectonic module
 URL: http://dune-project.org/
-Requires: dune-common dune-fufem dune-tnnmg
-#Libs: -L${libdir} -ldunetectonic
+Requires: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Libs: -L${libdir}
 Cflags: -I${includedir}

dune.module

@@ -2,9 +2,7 @@
 # Dune module information file #
 ################################
 
-#Name of the module
 Module: dune-tectonic
-Version: 2.1-1
-Maintainer: pipping@mi.fu-berlin.de
-#depending on 
-Depends: dune-common dune-fufem dune-tnnmg
+Version: 2.5-1
+Maintainer: elias.pipping@fu-berlin.de
+Depends: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid

dune/CMakeLists.txt

+add_subdirectory(tectonic)

dune/Makefile.am

-SUBDIRS = tectonic
-
-include $(top_srcdir)/am/global-rules

dune/tectonic/CMakeLists.txt

+install(FILES
+  body.hh
+  frictiondata.hh
+  frictionpotential.hh
+  globalfrictiondata.hh
+  globalfriction.hh
+  globalratestatefriction.hh
+  gravity.hh
+  localfriction.hh
+  minimisation.hh
+  myblockproblem.hh
+  mydirectionalconvexfunction.hh
+  quadraticenergy.hh
+  tectonic.hh
+  DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/Makefile.am

-tectonicincludedir = $(includedir)/dune/tectonic
-tectonicinclude_HEADERS = tectonic.hh
-
-include $(top_srcdir)/am/global-rules

dune/tectonic/body.hh

@@ -2,8 +2,9 @@
 #define DUNE_TECTONIC_BODY_HH
 
 template <int dimension> struct Body {
-  using ScalarFunction = Dune::VirtualFunction<
-      Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
+  using ScalarFunction =
+      Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                            Dune::FieldVector<double, 1>>;
   using VectorField =
       Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
                             Dune::FieldVector<double, dimension>>;

dune/tectonic/closestpointprojection.hh

-// Based on the class from dune-contact
-
-#ifndef CLOSEST_POINT_PROJECTION_HH
-#define CLOSEST_POINT_PROJECTION_HH
-
-template <class Coordinate> struct BoundarySegment {
-  int nVertices;
-  std::vector<Coordinate> vertices;
-};
-
-template <class Coordinate>
-Coordinate computeClosestPoint(const Coordinate &target,
-                               const BoundarySegment<Coordinate> &segment,
-                               int gaussSeidelIterations) {
-  /**
-   *  The functional to be minimized is: 0.5*norm(\sum_i lambda_i corner_i  -
-   *target)^2
-   *  where lambda_i are barycentric coordinates, i.e. 0\leq lambda_i \leq 1 and
-   *\sum_i lambda_i=1
-   *
-   *  The resulting quadratic minimization problem is given by 0.5
-   *Lambda^T*A*Lambda - l^T*Lambda
-   *  with
-   *      A_ij = (corner_i,corner_j)  and l_i = (corner_i, target)
-  */
-  int nCorners = segment.nVertices;
-  double A[nCorners][nCorners];
-
-  for (int i = 0; i < nCorners; i++)
-    for (int j = 0; j <= i; j++)
-      A[i][j] = segment.vertices[i] * segment.vertices[j];
-
-  // make symmetric
-  for (int i = 0; i < nCorners; i++)
-    for (int j = i + 1; j < nCorners; j++)
-      A[i][j] = A[j][i];
-
-  std::vector<double> l(nCorners);
-  for (int i = 0; i < nCorners; i++)
-    l[i] = target * segment.vertices[i];
-
-  // choose a feasible initial solution
-  std::vector<double> sol(nCorners);
-  for (int i = 0; i < nCorners; i++)
-    sol[i] = 1.0 / nCorners;
-
-  // use a Gauß-Seidel like method for possibly non-convex quadratic problems
-  // with simplex constraints.
-  for (int i = 0; i < gaussSeidelIterations; i++) {
-
-    // compute the residual
-    std::vector<double> rhs = l;
-
-    // compute rhs -= A*sol_
-    for (int k = 0; k < nCorners; k++)
-      for (int j = 0; j < nCorners; j++)
-        rhs[k] -= A[k][j] * sol[j];
-
-    // use the edge vectors as search directions
-    double alpha = 0.0;
-    for (int j1 = 0; j1 < (nCorners - 1); ++j1)
-      for (int j2 = j1 + 1; j2 < nCorners; ++j2) {
-        // compute matrix entry and rhs for edge direction
-        double a = A[j1][j1] - A[j1][j2] - A[j2][j1] + A[j2][j2];
-        double b = rhs[j1] - rhs[j2];
-
-        // compute minimizer for correction
-        if (a > 0) {
-          alpha = b / a;
-          // project alpha such that we stay positiv
-          if ((sol[j2] - alpha) < -1e-14) {
-            alpha = sol[j2];
-          } else if ((alpha < -1e-14) && ((sol[j1] + alpha) < -1e-14)) {
-            alpha = -sol[j1];
-          }
-        } else {
-          // if concave the minimum is achieved at one of the boundaries,
-          // i.e. at sol_[j2] or -sol_[j1]
-
-          double sum = sol[j1] + sol[j2];
-
-          double lValue = 0.5 * A[j1][j1] * sum * sum - rhs[j1] * sum;
-          double uValue = 0.5 * A[j2][j2] * sum * sum - rhs[j2] * sum;
-
-          alpha = (lValue < uValue) ? sol[j2] : -sol[j1];
-        }
-        // apply correction
-        sol[j1] += alpha;
-        sol[j2] -= alpha;
-
-        // update the local residual for corrections
-        for (int p = 0; p < nCorners; p++)
-          rhs[p] -= (A[p][j1] - A[p][j2]) * alpha;
-      }
-  }
-
-  double eps = 1e-5;
-  double sum = 0;
-  for (size_t i = 0; i < sol.size(); i++) {
-    if (sol[i] < -eps)
-      DUNE_THROW(Dune::Exception, "No barycentric coords " << sol[1] << " nr. "
-                                                           << i);
-    sum += sol[i];
-  }
-  if (sum > 1 + eps)
-    DUNE_THROW(Dune::Exception, "No barycentric coords, sum: " << sum);
-
-  Coordinate projected(0);
-  for (int i = 0; i < nCorners; i++)
-    projected.axpy(sol[i], segment.vertices[i]);
-
-  return projected;
-}
-#endif

dune/tectonic/frictionpotential.hh

@@ -29,35 +29,36 @@ class FrictionPotential {
 
 class TruncatedRateState : public FrictionPotential {
 public:
-  TruncatedRateState(double coefficient, double _normalStress, FrictionData _fd)
-      : fd(_fd), weight(coefficient), normalStress(_normalStress) {}
+  TruncatedRateState(double _weight, double _weightedNormalStress,
+                     FrictionData _fd)
+      : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
 
-  double coefficientOfFriction(double V) const {
+  double coefficientOfFriction(double V) const override {
     if (V <= Vmin)
       return 0.0;
 
     return fd.a * std::log(V / Vmin);
   }
 
-  double differential(double V) const {
-    return weight * (fd.C - normalStress * coefficientOfFriction(V));
+  double differential(double V) const override {
+    return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
   }
 
-  double second_deriv(double V) const {
+  double second_deriv(double V) const override {
     if (V <= Vmin)
       return 0;
 
-    return weight * (-normalStress) * (fd.a / V);
+    return -weightedNormalStress * (fd.a / V);
   }
 
-  double regularity(double V) const {
+  double regularity(double V) const override {
     if (std::abs(V - Vmin) < 1e-14) // TODO
       return std::numeric_limits<double>::infinity();
 
     return std::abs(second_deriv(V));
   }
 
-  void updateAlpha(double alpha) {
+  void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
   }
@@ -65,31 +66,33 @@ class TruncatedRateState : public FrictionPotential {
 private:
   FrictionData const fd;
   double const weight;
-  double const normalStress;
+  double const weightedNormalStress;
   double Vmin;
 };
 
 class RegularisedRateState : public FrictionPotential {
 public:
-  RegularisedRateState(double coefficient, double _normalStress,
+  RegularisedRateState(double _weight, double _weightedNormalStress,
                        FrictionData _fd)
-      : fd(_fd), weight(coefficient), normalStress(_normalStress) {}
+      : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
 
-  double coefficientOfFriction(double V) const {
+  double coefficientOfFriction(double V) const override {
     return fd.a * std::asinh(0.5 * V / Vmin);
   }
 
-  double differential(double V) const {
-    return weight * (fd.C - normalStress * coefficientOfFriction(V));
+  double differential(double V) const override {
+    return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
   }
 
-  double second_deriv(double V) const {
-    return weight * (-normalStress) * fd.a / std::hypot(2.0 * Vmin, V);
+  double second_deriv(double V) const override {
+    return -weightedNormalStress * fd.a / std::hypot(2.0 * Vmin, V);
   }
 
-  double regularity(double V) const { return std::abs(second_deriv(V)); }
+  double regularity(double V) const override {
+    return std::abs(second_deriv(V));
+  }
 
-  void updateAlpha(double alpha) {
+  void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
   }
@@ -97,22 +100,22 @@ class RegularisedRateState : public FrictionPotential {
 private:
   FrictionData const fd;
   double const weight;
-  double const normalStress;
+  double const weightedNormalStress;
   double Vmin;
 };
 
-class TrivialFunction : public FrictionPotential {
+class ZeroFunction : public FrictionPotential {
 public:
-  double evaluate(double) const { return 0; }
+  double evaluate(double) const override { return 0; }
 
-  double coefficientOfFriction(double s) const { return 0; }
+  double coefficientOfFriction(double s) const override { return 0; }
 
-  double differential(double) const { return 0; }
+  double differential(double) const override { return 0; }
 
-  double second_deriv(double) const { return 0; }
+  double second_deriv(double) const override { return 0; }
 
-  double regularity(double) const { return 0; }
+  double regularity(double) const override { return 0; }
 
-  void updateAlpha(double) {}
+  void updateAlpha(double) override {}
 };
 #endif

dune/tectonic/geocoordinate.hh

+#ifndef DUNE_TECTONIC_GEOCOORDINATE_HH
+#define DUNE_TECTONIC_GEOCOORDINATE_HH
+
+// tiny helper to make a common piece of code pleasanter to read
+
+template <class Geometry>
+typename Geometry::GlobalCoordinate geoToPoint(Geometry geo) {
+  assert(geo.corners() == 1);
+  return geo.corner(0);
+}
+
+#endif