diff --git a/.gitignore b/.gitignore
index c68e2924dce682170ab15fa77236298ad6514549..f47390d21927c4cf3aa0ca9a910d12fc7c6ea7a1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -26,4 +26,5 @@
/test-driver
Makefile
Makefile.in
+src/sand-wedge
src/sliding-block-?D
diff --git a/configure.ac b/configure.ac
index 0ff34a5477e8725a892f83f2efad9633f779b1ff..f7b98bc035b8f8ddaee0e3896b24fc72f4278c7d 100644
--- a/configure.ac
+++ b/configure.ac
@@ -9,6 +9,16 @@ 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
diff --git a/m4/dune-tectonic.m4 b/m4/dune-tectonic.m4
index 9bde6699ff6e14cd1321ec36846d9a1fd607bfc9..c8a3afe84ef8759a4c10a7f4e8177076cf89ace2 100644
--- a/m4/dune-tectonic.m4
+++ b/m4/dune-tectonic.m4
@@ -6,7 +6,9 @@ dnl -*- autoconf -*-
# Additional checks needed to build dune-tectonic
# This macro should be invoked by every module which depends on dune-tectonic, as
# well as by dune-tectonic itself
-AC_DEFUN([DUNE_TECTONIC_CHECKS])
+AC_DEFUN([DUNE_TECTONIC_CHECKS],[
+ AC_REQUIRE([AX_BOOST_BASE])
+])
# Additional checks needed to find dune-tectonic
# This macro should be invoked by every module which depends on dune-tectonic, but
diff --git a/src/Makefile.am b/src/Makefile.am
index b8a7ad36909d3695602b0aeb46dbc0b66353d2f8..91d1688e3f719530a01602347bd5ddabbb8e8bc8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS = sliding-block-2D
+bin_PROGRAMS = sand-wedge sliding-block-2D
common_sources = \
assemblers.cc \
@@ -8,6 +8,10 @@ common_sources = \
timestepping.cc \
vtk.cc
+sand_wedge_SOURCES = $(common_sources) sand-wedge.cc
+sand_wedge_CPPFLAGS = \
+ $(AM_CPPFLAGS) $(BOOST_CPPFLAGS) \
+ -Ddatadir=\"$(abs_srcdir)/sand-wedge-data/\" -DDIM=2
sliding_block_2D_SOURCES = $(common_sources) sliding-block.cc
sliding_block_2D_CPPFLAGS = \
$(AM_CPPFLAGS) \
@@ -41,4 +45,9 @@ AM_LDFLAGS = \
$(ALUGRID_LDFLAGS) \
$(PYTHON_LDFLAGS)
+if CAIROMM
+AM_CPPFLAGS += $(CAIROMM_CFLAGS)
+LDADD += $(CAIROMM_LIBS)
+endif
+
include $(top_srcdir)/am/global-rules
diff --git a/src/sand-wedge-data/boundaryconditions.py b/src/sand-wedge-data/boundaryconditions.py
new file mode 100644
index 0000000000000000000000000000000000000000..5143bd93f2d13a3a20fb02141756c13c7e2af32a
--- /dev/null
+++ b/src/sand-wedge-data/boundaryconditions.py
@@ -0,0 +1,16 @@
+class neumannCondition:
+ def __call__(self, relativeTime):
+ return 0
+
+class velocityDirichletCondition:
+ def __call__(self, relativeTime):
+ if relativeTime <= 0.25:
+ factor = 4*relativeTime;
+ else:
+ factor = 1
+ return -5e-5 * factor
+
+Functions = {
+ 'neumannCondition' : neumannCondition(),
+ 'velocityDirichletCondition' : velocityDirichletCondition()
+}
diff --git a/src/sand-wedge-data/geometry.tex b/src/sand-wedge-data/geometry.tex
new file mode 100644
index 0000000000000000000000000000000000000000..7b2ca811766ca5756771d7215cb73dca3652d68a
--- /dev/null
+++ b/src/sand-wedge-data/geometry.tex
@@ -0,0 +1,77 @@
+\documentclass[tikz]{minimal}
+
+\usepackage{tikz}
+\usetikzlibrary{calc}
+\usetikzlibrary{decorations.pathreplacing}
+
+
+\usepackage{siunitx}
+
+\begin{document}
+\begin{tikzpicture}[scale=12]
+ \pgfmathsetmacro{\hypot}{sqrt(1+.27^2)};
+ \pgfmathsetmacro{\ssin}{sin(15)};
+ \pgfmathsetmacro{\scos}{cos(15)};
+ \pgfmathsetmacro{\lsin}{sin(75)};
+
+ \coordinate (A) at (0,0);
+ \node at (A) [left] {A};
+ \coordinate (B) at (\hypot,0);
+ \node at (B) [right] {B};
+ \coordinate (C) at (\lsin, .27*\lsin);
+ \node at (C) [above] {C};
+
+ \draw (A)
+ -- (B)
+ -- (C)
+ -- node [above=.5cm, sloped] {$\overline{AC}=\SI{100}{cm}$} (A);
+
+ \coordinate (Y) at ($(0.35/\lsin, 0)$);
+ \node at (Y) [below right] {Y};
+ \coordinate (Z) at ($(0.35*\lsin,0.35*\ssin)$);
+ \node at (Z) [above] {Z};
+ \coordinate (X) at ($(Y) - (0.20,0)$);
+ \node at (X) [below left] {X};
+ \path let \p1 = (X), \p2 = (Y), \p3 = (Z) in
+ coordinate (U) at ($(\x1 * \x3 / \x2, \x1 * \y3 / \x2)$);
+ \node at (U) [above] {U};
+
+ \path (A) -- node [above=.25cm, sloped] {$\overline{AZ} = \SI{35}{cm}$} (Z);
+
+ \draw[color=red, thick] (X)
+ -- node [below=.25cm, sloped] {$\overline{XY}=\SI{20}{cm}$} (Y);
+ \draw[dashed] (Y) -- (Z);
+ \draw[dashed] (U) -- (X);
+
+ \coordinate (K) at ($(B) - (0.06 / \ssin,0)$);
+ \node at (K) [below] {K};
+ \coordinate (M) at ($(B) + (-0.06 * \ssin, 0.06 * \lsin)$);
+ \node at (M) [right] {M};
+
+ \path (C) -- node [right=.5cm] {$\overline{CM} = \SI{21}{cm}$} (M);
+
+ \draw[thick] (K)
+ -- (B)
+ -- node [right=.75cm] {$\overline{BM}=\SI{6}{cm}$} (M) -- cycle;
+
+ \coordinate (G) at ($(A) ! 0.5 ! (X)$);
+ \node at (G) [below] {G};
+ \coordinate (H) at ($(X) ! 0.5 ! (Y)$);
+ \node at (H) [below] {H};
+ \coordinate (J) at ($(Y) ! 0.5 ! (B)$);
+ \node at (J) [below] {J};
+
+ \coordinate (I) at ($(Y) + (G)$);
+ \node at (I) [below] {I};
+
+ % Note: Lengths and position are arbitrarily chosen here
+ \draw[->] (0.4,0.45) -- node [above, sloped] {vertical} +(-0.15*\ssin,0.15*\scos);
+ \draw[->] (0.4,0.45) -- node [above, sloped] {horizontal} +(-0.15*\scos,-0.15*\ssin);
+
+ \node[align=left] at (0.5,-0.125) {
+ $Z$: coast line\\
+ $\overline{XY}$: velocity weakening zone\\
+ $BKM$: visco-elastic domain};
+\end{tikzpicture}
+
+\end{document}
diff --git a/src/sand-wedge-data/mybody.hh b/src/sand-wedge-data/mybody.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a40f79a6cc8436b291c6d5a8e8c7214423cb0a61
--- /dev/null
+++ b/src/sand-wedge-data/mybody.hh
@@ -0,0 +1,55 @@
+#ifndef SRC_SAND_WEDGE_DATA_MYBODY_HH
+#define SRC_SAND_WEDGE_DATA_MYBODY_HH
+
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/functions/constantfunction.hh>
+
+#include <dune/tectonic/body.hh>
+#include <dune/tectonic/gravity.hh>
+
+#include "twopiece.hh"
+#include "mygeometry.hh"
+
+template <int dimension> class MyBody : public Body<dimension> {
+ using typename Body<dimension>::ScalarFunction;
+ using typename Body<dimension>::VectorField;
+
+public:
+ MyBody(Dune::ParameterTree const &parset)
+ : poissonRatio_(parset.get<double>("body.poissonRatio")),
+ youngModulus_(3.0 * parset.get<double>("body.bulkModulus") *
+ (1.0 - 2.0 * poissonRatio_)),
+ shearViscosityField_(
+ parset.get<double>("body.elastic.shearViscosity"),
+ parset.get<double>("body.viscoelastic.shearViscosity")),
+ bulkViscosityField_(
+ parset.get<double>("body.elastic.bulkViscosity"),
+ parset.get<double>("body.viscoelastic.bulkViscosity")),
+ densityField_(parset.get<double>("body.elastic.density"),
+ parset.get<double>("body.viscoelastic.density")),
+ gravityField_(densityField_, MyGeometry::zenith,
+ parset.get<double>("gravity")) {}
+
+ double getPoissonRatio() const override { return poissonRatio_; }
+ double getYoungModulus() const override { return youngModulus_; }
+ ScalarFunction const &getShearViscosityField() const override {
+ return shearViscosityField_;
+ }
+ ScalarFunction const &getBulkViscosityField() const override {
+ return bulkViscosityField_;
+ }
+ ScalarFunction const &getDensityField() const override {
+ return densityField_;
+ }
+ VectorField const &getGravityField() const override { return gravityField_; }
+
+private:
+ double const poissonRatio_;
+ double const youngModulus_;
+ TwoPieceFunction const shearViscosityField_;
+ TwoPieceFunction const bulkViscosityField_;
+ TwoPieceFunction const densityField_;
+ Gravity<dimension> const gravityField_;
+};
+#endif
diff --git a/src/sand-wedge-data/mygeometry.cc b/src/sand-wedge-data/mygeometry.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5629772f2299b5e2fab4bb7449541e98d31cf0ec
--- /dev/null
+++ b/src/sand-wedge-data/mygeometry.cc
@@ -0,0 +1,161 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <fstream>
+
+#include "mygeometry.hh"
+
+double MyGeometry::verticalProjection(LocalVector const &x) {
+ return zenith * x;
+}
+double MyGeometry::horizontalProjection(LocalVector const &x) {
+ return rotatedZenith * x;
+}
+
+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;
+}
+
+void MyGeometry::render() {
+#ifdef HAVE_CAIROMM
+ std::string const filename = "geometry.png";
+ double const width = 600;
+ double const height = 400;
+ double const widthScale = 400;
+ double const heightScale = 400;
+
+ auto surface =
+ Cairo::ImageSurface::create(Cairo::FORMAT_ARGB32, width, height);
+ auto cr = Cairo::Context::create(surface);
+
+ auto const setRGBColor = [&](int colour) {
+ cr->set_source_rgb(((colour & 0xFF0000) >> 16) / 255.0,
+ ((colour & 0x00FF00) >> 8) / 255.0,
+ ((colour & 0x0000FF) >> 0) / 255.0);
+ };
+ auto const moveTo = [&](LocalVector const &v) { cr->move_to(v[0], -v[1]); };
+ auto const lineTo = [&](LocalVector const &v) { cr->line_to(v[0], -v[1]); };
+
+ cr->scale(widthScale, heightScale);
+ cr->translate(0.1, 0.1);
+ cr->rotate(leftAngle);
+ cr->set_line_width(0.0025);
+
+ // triangle
+ {
+ moveTo(A);
+ lineTo(B);
+ lineTo(C);
+ cr->close_path();
+ cr->stroke();
+ }
+
+ // dashed lines
+ {
+ cr->save();
+ std::vector<double> dashPattern = { 0.005 };
+ cr->set_dash(dashPattern, 0);
+ moveTo(Z);
+ lineTo(Y);
+ moveTo(U);
+ lineTo(X);
+ cr->stroke();
+ cr->restore();
+ }
+
+ // fill viscoelastic region
+ {
+ cr->save();
+ setRGBColor(0x0097E0);
+ moveTo(B);
+ lineTo(K);
+ lineTo(M);
+ cr->fill();
+ cr->restore();
+ }
+
+ // mark weakening region
+ {
+ cr->save();
+ setRGBColor(0x7AD3FF);
+ cr->set_line_width(0.005);
+ moveTo(X);
+ lineTo(Y);
+ cr->stroke();
+ cr->restore();
+ }
+
+ // mark points
+ {
+ auto const drawCircle = [&](LocalVector const &v) {
+ cr->arc(v[0], -v[1], 0.0075, -M_PI, M_PI); // x,y,radius,angle1,angle2
+ cr->fill();
+ };
+
+ cr->save();
+ setRGBColor(0x002F47);
+ drawCircle(A);
+ drawCircle(B);
+ drawCircle(C);
+ drawCircle(Y);
+ drawCircle(X);
+ drawCircle(Z);
+ drawCircle(U);
+ drawCircle(K);
+ drawCircle(M);
+ drawCircle(G);
+ drawCircle(H);
+ drawCircle(J);
+ drawCircle(I);
+ cr->restore();
+ }
+
+ // labels
+ {
+ auto const label = [&](LocalVector const &v,
+ std::string l) { // TODO: implicit?
+ moveTo(v);
+ cr->rel_move_to(0.005, -0.02);
+ cr->show_text(l);
+ };
+ auto font = Cairo::ToyFontFace::create(
+ "monospace", Cairo::FONT_SLANT_NORMAL, Cairo::FONT_WEIGHT_NORMAL);
+
+ cr->save();
+ cr->set_font_face(font);
+ cr->set_font_size(0.03);
+
+ label(A, "A");
+ label(B, "B");
+ label(C, "C");
+ label(K, "K");
+ label(M, "M");
+ label(U, "U");
+ label(X, "X");
+ label(Y, "Y");
+ label(Z, "Z");
+ label(G, "G");
+ label(H, "H");
+ label(J, "J");
+ label(I, "I");
+ cr->restore();
+ }
+
+ surface->write_to_png(filename);
+#endif
+}
diff --git a/src/sand-wedge-data/mygeometry.hh b/src/sand-wedge-data/mygeometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d9955b2c4d0c5039306da116eeaf0aa50ce1224a
--- /dev/null
+++ b/src/sand-wedge-data/mygeometry.hh
@@ -0,0 +1,82 @@
+#ifndef MY_GEOMETRY_HH
+#define MY_GEOMETRY_HH
+
+#include <dune/common/fassign.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/arithmetic.hh>
+
+#ifdef HAVE_CAIROMM
+#include <cairomm/context.h>
+#include <cairomm/fontface.h>
+#include <cairomm/surface.h>
+#endif
+
+namespace MyGeometry {
+namespace {
+ using LocalVector = Dune::FieldVector<double, 2>;
+
+ // kludge because fieldvectors have no initialiser_list constructor, see
+ // https://dune-project.org/flyspray/index.php?do=details&task_id=1166
+ LocalVector generateVector(double x, double y) {
+ LocalVector tmp;
+ tmp <<= x, y;
+ return tmp;
+ }
+
+ double degreesToRadians(double degrees) {
+ return M_PI * degrees / 180;
+ };
+
+ LocalVector convexCombination(LocalVector const &x, LocalVector const &y) {
+ LocalVector ret;
+ ret = 0.0;
+ Arithmetic::addProduct(ret, 0.5, x);
+ Arithmetic::addProduct(ret, 0.5, y);
+ return ret;
+ }
+
+ double leftLeg = 1.00;
+ double rightLeg = 0.27;
+ double leftAngle = degreesToRadians(15);
+ double rightAngle = degreesToRadians(75);
+
+ double const depth = 0.10;
+}
+
+LocalVector const A = generateVector(0, 0);
+LocalVector const B = generateVector(std::hypot(leftLeg, rightLeg), 0);
+LocalVector const C = generateVector(leftLeg * std::sin(rightAngle),
+ leftLeg * std::sin(leftAngle));
+
+LocalVector const Y = generateVector(0.35 / std::sin(rightAngle), 0);
+LocalVector const X = generateVector(Y[0] - 0.20, 0);
+LocalVector const Z =
+ generateVector(0.35 * std::sin(rightAngle), 0.35 * std::sin(leftAngle));
+
+LocalVector const U = generateVector(Z[0] * X[0] / Y[0], Z[1] * X[0] / Y[0]);
+
+LocalVector const K = generateVector(B[0] - 0.06 / std::sin(leftAngle), 0);
+LocalVector const M = generateVector(B[0] - 0.06 * std::sin(leftAngle),
+ 0.06 * std::sin(rightAngle));
+
+LocalVector const G = convexCombination(A, X);
+LocalVector const H = convexCombination(X, Y);
+LocalVector const J = convexCombination(Y, B);
+
+LocalVector const I = generateVector(Y[0] + G[0], 0);
+
+LocalVector const zenith =
+ generateVector(-std::sin(leftAngle), std::cos(leftAngle));
+
+LocalVector const rotatedZenith =
+ generateVector(-std::cos(leftAngle), -std::sin(leftAngle));
+
+double verticalProjection(LocalVector const &);
+double horizontalProjection(LocalVector const &);
+
+void write();
+
+void render();
+}
+#endif
diff --git a/src/sand-wedge-data/myglobalfrictiondata.hh b/src/sand-wedge-data/myglobalfrictiondata.hh
new file mode 100644
index 0000000000000000000000000000000000000000..64db619cb3d6dda72fa80c85070ffa6ef50dfd62
--- /dev/null
+++ b/src/sand-wedge-data/myglobalfrictiondata.hh
@@ -0,0 +1,43 @@
+#ifndef SRC_SAND_WEDGE_DATA_MYGLOBALFRICTIONDATA_HH
+#define SRC_SAND_WEDGE_DATA_MYGLOBALFRICTIONDATA_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "mygeometry.hh"
+#include "patchfunction.hh"
+
+template <int dimension>
+class MyGlobalFrictionData : public GlobalFrictionData<dimension> {
+private:
+ using typename GlobalFrictionData<dimension>::VirtualFunction;
+
+public:
+ MyGlobalFrictionData(Dune::ParameterTree const &parset)
+ : C_(parset.get<double>("C")),
+ L_(parset.get<double>("L")),
+ V0_(parset.get<double>("V0")),
+ a_(parset.get<double>("strengthening.a"),
+ parset.get<double>("weakening.a")),
+ b_(parset.get<double>("strengthening.b"),
+ parset.get<double>("weakening.b")),
+ mu0_(parset.get<double>("mu0")) {}
+
+ double virtual const &C() const override { return C_; }
+ double virtual const &L() const override { return L_; }
+ double virtual const &V0() const override { return V0_; }
+ VirtualFunction virtual const &a() const override { return a_; }
+ VirtualFunction virtual const &b() const override { return b_; }
+ double virtual const &mu0() const override { return mu0_; }
+
+private:
+ double const C_;
+ double const L_;
+ double const V0_;
+ PatchFunction const a_;
+ PatchFunction const b_;
+ double const mu0_;
+};
+#endif
diff --git a/src/sand-wedge-data/mygrid.cc b/src/sand-wedge-data/mygrid.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ec2343950395853711896cccd3c9a5183d23acbf
--- /dev/null
+++ b/src/sand-wedge-data/mygrid.cc
@@ -0,0 +1,76 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "mygrid.hh"
+
+template <class Grid> std::shared_ptr<Grid> constructGrid() {
+ Dune::GridFactory<Grid> gridFactory;
+
+ Dune::FieldMatrix<double, 3, DIM> vertices;
+ vertices[0] = MyGeometry::A;
+ vertices[1] = MyGeometry::B;
+ vertices[2] = MyGeometry::C;
+ for (size_t i = 0; i < vertices.N(); ++i)
+ gridFactory.insertVertex(vertices[i]);
+
+ Dune::GeometryType cell;
+ cell.makeTriangle();
+
+ std::vector<std::vector<unsigned int>> simplices = { { 0, 1, 2 } };
+
+ // sanity-check choices of simplices
+ for (size_t i = 0; i < simplices.size(); ++i) {
+ Dune::FieldMatrix<double, DIM, DIM> check;
+ for (size_t j = 0; j < DIM; ++j)
+ check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
+ assert(check.determinant() > 0);
+ gridFactory.insertElement(cell, simplices[i]);
+ }
+ return std::shared_ptr<Grid>(gridFactory.createGrid());
+}
+
+template <class GridView>
+MyFaces<GridView>::MyFaces(GridView const &gridView)
+ : lower(gridView), right(gridView) {
+ using Grid = typename GridView::Grid;
+ using LevelGridView = typename Grid::LevelGridView;
+ LevelGridView const rootView = gridView.grid().levelView(0);
+
+ lower.insertFacesByProperty([&](
+ typename Grid::LeafGridView::Intersection const &in) {
+ return isClose(0, in.geometry().center()[1]);
+ });
+
+ BoundaryPatch<LevelGridView> upperRoot(rootView);
+ upperRoot.insertFacesByProperty([&](
+ typename LevelGridView::Intersection const &in) {
+ auto const geometry = in.geometry();
+ auto const range = boost::irange(0, geometry.corners());
+ return boost::algorithm::all_of(range, [&](int i) {
+ auto const &co = geometry.corner(i);
+ return (isClose(co[1], MyGeometry::A[1]) &&
+ isClose(co[0], MyGeometry::A[0])) ||
+ (isClose(co[1], MyGeometry::C[1]) &&
+ isClose(co[0], MyGeometry::C[0]));
+ });
+ });
+ BoundaryPatchProlongator<Grid>::prolong(upperRoot, upper);
+
+ BoundaryPatch<LevelGridView> rightRoot(rootView);
+ rightRoot.insertFacesByProperty([&](
+ typename LevelGridView::Intersection const &in) {
+ auto const geometry = in.geometry();
+ auto const range = boost::irange(0, geometry.corners());
+ return boost::algorithm::all_of(range, [&](int i) {
+ auto const &co = geometry.corner(i);
+ return (isClose(co[1], MyGeometry::C[1]) &&
+ isClose(co[0], MyGeometry::C[0])) ||
+ (isClose(co[1], MyGeometry::B[1]) &&
+ isClose(co[0], MyGeometry::B[0]));
+ });
+ });
+ BoundaryPatchProlongator<Grid>::prolong(rightRoot, right);
+}
+
+#include "mygrid_tmpl.cc"
diff --git a/src/sand-wedge-data/mygrid.hh b/src/sand-wedge-data/mygrid.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1f02943dd6d3d415cbbb3cc3fe75419beb58b54a
--- /dev/null
+++ b/src/sand-wedge-data/mygrid.hh
@@ -0,0 +1,31 @@
+#ifndef SRC_SAND_WEDGE_DATA_MYGRID_HH
+#define SRC_SAND_WEDGE_DATA_MYGRID_HH
+
+#include <boost/range/irange.hpp>
+#include <boost/algorithm/cxx11/all_of.hpp>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wdeprecated-declarations"
+#include <dune/fufem/boundarypatch.hh>
+#pragma clang diagnostic pop
+#include <dune/fufem/boundarypatchprolongator.hh>
+
+#include "mygeometry.hh"
+
+template <class Grid> std::shared_ptr<Grid> constructGrid();
+
+template <class GridView> struct MyFaces {
+ BoundaryPatch<GridView> lower;
+ BoundaryPatch<GridView> right;
+ BoundaryPatch<GridView> upper;
+ MyFaces(GridView const &gridView);
+
+private:
+ bool isClose(double a, double b) {
+ return std::abs(a - b) < 1e-14;
+ };
+};
+#endif
diff --git a/src/sand-wedge-data/mygrid_tmpl.cc b/src/sand-wedge-data/mygrid_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..3afed427fc51b0dd542d6a01e316fc4b5ce35d72
--- /dev/null
+++ b/src/sand-wedge-data/mygrid_tmpl.cc
@@ -0,0 +1,9 @@
+#ifndef DIM
+#error DIM unset
+#endif
+
+#include "explicitgrid.hh"
+
+template std::shared_ptr<Grid> constructGrid();
+
+template struct MyFaces<GridView>;
diff --git a/src/sand-wedge-data/parset.cfg b/src/sand-wedge-data/parset.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..4287864a84c0f93bf481637cc87412d3996aa13f
--- /dev/null
+++ b/src/sand-wedge-data/parset.cfg
@@ -0,0 +1,78 @@
+# -*- mode:conf -*-
+gravity = 9.81 # [m/s^2]
+
+[io]
+printProgress = false
+writeVTK = false
+
+[problem]
+finalTime = 1800 # [s]
+
+[body]
+bulkModulus = 1e5 # [Pa]
+poissonRatio = 0.3 # [1] 0.2 - 0.3
+[body.elastic]
+density = 900 # [kg/m^3]
+shearViscosity = 1e3 # [Pas] 0
+bulkViscosity = 1e3 # [Pas] 0
+[body.viscoelastic]
+density = 1000 # [kg/m^3]
+shearViscosity = 1e4 # [Pas]
+bulkViscosity = 1e4 # [Pas]
+
+[boundary.friction]
+C = 10 # [Pa]
+mu0 = 0.7 # [1]
+V0 = 5e-5 # [m/s]
+L = 2e-5 # [m] ?
+initialLogState = 0 # [ ] ?
+stateModel = Dieterich
+[boundary.friction.weakening]
+a = 0.015 # [1] ?
+b = 0.030 # [1] ?
+[boundary.friction.strengthening]
+a = 0.030 # [1] ?
+b = 0.015 # [1] ?
+
+[timeSteps]
+number = 100000
+scheme = newmark
+
+[grid]
+refinements = 4
+
+[u0.solver]
+tolerance = 1e-10
+maximumIterations = 100000
+verbosity = quiet
+
+[a0.solver]
+tolerance = 1e-10
+maximumIterations = 100000
+verbosity = quiet
+
+[v.solver]
+tolerance = 1e-10
+maximumIterations = 100000
+verbosity = quiet
+
+[v.fpi]
+tolerance = 1e-10
+maximumIterations = 10000
+relaxation = 0.5
+requiredReduction = 0.5
+
+[solver.tnnmg.linear]
+maxiumumIterations = 100000
+tolerance = 1e-10
+pre = 3
+cycle = 1 # 1 = V, 2 = W, etc.
+post = 3
+
+[solver.tnnmg.main]
+pre = 1
+multi = 5 # number of multigrid steps
+post = 0
+
+[localsolver]
+steps = 100
diff --git a/src/sand-wedge-data/patchfunction.hh b/src/sand-wedge-data/patchfunction.hh
new file mode 100644
index 0000000000000000000000000000000000000000..871e65fdb1fc653261d6846d28a39e79cea30b70
--- /dev/null
+++ b/src/sand-wedge-data/patchfunction.hh
@@ -0,0 +1,40 @@
+#ifndef SRC_SAND_WEDGE_DATA_PATCHFUNCTION_HH
+#define SRC_SAND_WEDGE_DATA_PATCHFUNCTION_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+class PatchFunction
+ : public Dune::VirtualFunction<Dune::FieldVector<double, 2>,
+ Dune::FieldVector<double, 1>> {
+private:
+ bool static isBetween(double x, double lower, double upper) {
+ return lower <= x and x <= upper;
+ }
+
+ bool static isClose(double a, double b) {
+ return std::abs(a - b) < 1e-14;
+ };
+
+ bool insideRegion2(Dune::FieldVector<double, 2> const &z) const {
+ assert(isClose(0, z[1]));
+ return isBetween(z[0], _X[0], _Y[0]);
+ }
+
+ Dune::FieldVector<double, 2> const &_X;
+ Dune::FieldVector<double, 2> const &_Y;
+
+ double const _v1;
+ double const _v2;
+
+public:
+ PatchFunction(double v1, double v2)
+ : _X(MyGeometry::X), _Y(MyGeometry::Y), _v1(v1), _v2(v2) {}
+
+ void evaluate(Dune::FieldVector<double, 2> const &x,
+ Dune::FieldVector<double, 1> &y) const {
+ y = insideRegion2(x) ? _v2 : _v1;
+ }
+};
+#endif
diff --git a/src/sand-wedge-data/twopiece.hh b/src/sand-wedge-data/twopiece.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f243aa0fb10aaec8e68b935f12b25d6b4103310d
--- /dev/null
+++ b/src/sand-wedge-data/twopiece.hh
@@ -0,0 +1,38 @@
+#ifndef SRC_SAND_WEDGE_DATA_TWOPIECE_HH
+#define SRC_SAND_WEDGE_DATA_TWOPIECE_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include "mygeometry.hh"
+
+class TwoPieceFunction
+ : public Dune::VirtualFunction<Dune::FieldVector<double, 2>,
+ Dune::FieldVector<double, 1>> {
+private:
+ bool liesBelow(Dune::FieldVector<double, 2> const &x,
+ Dune::FieldVector<double, 2> const &y,
+ Dune::FieldVector<double, 2> const &z) const {
+ return (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1] - x[1];
+ };
+ bool insideRegion2(Dune::FieldVector<double, 2> const &z) const {
+ return liesBelow(_K, _M, z);
+ };
+
+ Dune::FieldVector<double, 2> const &_K;
+ Dune::FieldVector<double, 2> const &_M;
+
+ double const _v1;
+ double const _v2;
+
+public:
+ TwoPieceFunction(double v1, double v2)
+ : _K(MyGeometry::K), _M(MyGeometry::M), _v1(v1), _v2(v2) {}
+
+ void evaluate(Dune::FieldVector<double, 2> const &x,
+ Dune::FieldVector<double, 1> &y) const {
+ y = insideRegion2(x) ? _v2 : _v1;
+ }
+};
+#endif
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e8ae93c683fb840477c3f66b2393383967c404d3
--- /dev/null
+++ b/src/sand-wedge.cc
@@ -0,0 +1,544 @@
+#include <Python.h>
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifndef HAVE_PYTHON
+#error Python is required
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#ifndef datadir
+#error datadir unset
+#endif
+
+#ifndef DIM
+#error DIM unset
+#endif
+
+#if !HAVE_ALUGRID
+#error ALUGRID is required
+#endif
+
+#include <cmath>
+#include <exception>
+#include <fstream>
+#include <iostream>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wignored-qualifiers"
+#include <dune/grid/alugrid.hh>
+#pragma clang diagnostic pop
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wdeprecated-declarations"
+#include <dune/fufem/boundarypatch.hh>
+#pragma clang diagnostic pop
+
+#include <dune/fufem/dunepython.hh>
+#include <dune/fufem/functions/basisgridfunction.hh>
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/sharedpointermap.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/norms/sumnorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/solvers/solver.hh>
+#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
+#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
+#include <dune/tnnmg/problem-classes/convexproblem.hh>
+
+#include <dune/tectonic/myblockproblem.hh>
+#include <dune/tectonic/globalnonlinearity.hh>
+
+#include "assemblers.hh"
+#include "enum_parser.cc"
+#include "enum_scheme.cc"
+#include "enum_state_model.cc"
+#include "enum_verbosity.cc"
+#include "enums.hh"
+#include "friction_writer.hh"
+#include "sand-wedge-data/mybody.hh"
+#include "sand-wedge-data/mygeometry.hh"
+#include "sand-wedge-data/mygeometry.cc" // FIXME
+#include "sand-wedge-data/myglobalfrictiondata.hh"
+#include "sand-wedge-data/mygrid.cc" // FIXME
+#include "sand-wedge-data/mygrid.hh"
+#include "solverfactory.hh"
+#include "state.hh"
+#include "timestepping.hh"
+#include "vtk.hh"
+
+size_t const dims = DIM;
+
+void initPython() {
+ Python::start();
+
+ Python::run("import sys");
+ Python::run("sys.path.append('" datadir "')");
+}
+
+template <class GridView> class SpecialWriter {
+ using LocalVector = Dune::FieldVector<double, dims>;
+ using Element = typename GridView::Grid::template Codim<0>::Entity;
+ using ElementPointer =
+ typename GridView::Grid::template Codim<0>::EntityPointer;
+
+ void writeMagnitude(LocalVector const &v) {
+ writer_ << v[0] << " "; // Because of Dirichlet conditions, this is like a
+ // signed norm
+ }
+
+ void writeHorizontal(LocalVector const &v) {
+ writer_ << MyGeometry::horizontalProjection(v) << " ";
+ }
+ void writeVertical(LocalVector const &v) {
+ writer_ << MyGeometry::verticalProjection(v) << " ";
+ }
+
+ std::pair<ElementPointer, LocalVector> globalToLocal(LocalVector const &x)
+ const {
+ auto const element = hsearch_.findEntity(x);
+ return std::pair<ElementPointer, LocalVector>(element,
+ element->geometry().local(x));
+ }
+
+ std::fstream writer_;
+
+ Dune::HierarchicSearch<typename GridView::Grid, GridView> const hsearch_;
+
+ std::pair<ElementPointer, LocalVector> const G;
+ std::pair<ElementPointer, LocalVector> const H;
+ std::pair<ElementPointer, LocalVector> const J;
+ std::pair<ElementPointer, LocalVector> const I;
+ std::pair<ElementPointer, LocalVector> const U;
+ std::pair<ElementPointer, LocalVector> const Z;
+
+public:
+ SpecialWriter(std::string filename, GridView const &gridView)
+ : writer_(filename, std::fstream::out),
+ hsearch_(gridView.grid(), gridView),
+ G(globalToLocal(MyGeometry::G)),
+ H(globalToLocal(MyGeometry::H)),
+ J(globalToLocal(MyGeometry::J)),
+ I(globalToLocal(MyGeometry::I)),
+ U(globalToLocal(MyGeometry::U)),
+ Z(globalToLocal(MyGeometry::Z)) {
+ writer_ << "|G| |H| |J| |I| Uv Uh Zv Zh" << std::endl;
+ }
+
+ void write(VirtualGridFunction<typename GridView::Grid, LocalVector> const &
+ specialField) {
+ LocalVector value;
+
+ specialField.evaluateLocal(*G.first, G.second, value);
+ writeMagnitude(value);
+
+ specialField.evaluateLocal(*H.first, H.second, value);
+ writeMagnitude(value);
+
+ specialField.evaluateLocal(*J.first, J.second, value);
+ writeMagnitude(value);
+
+ specialField.evaluateLocal(*I.first, I.second, value);
+ writeMagnitude(value);
+
+ specialField.evaluateLocal(*U.first, U.second, value);
+ writeVertical(value);
+ writeHorizontal(value);
+
+ specialField.evaluateLocal(*Z.first, Z.second, value);
+ writeVertical(value);
+ writeHorizontal(value);
+
+ writer_ << std::endl;
+ }
+};
+
+int main(int argc, char *argv[]) {
+ try {
+ Dune::ParameterTree parset;
+ Dune::ParameterTreeParser::readINITree(datadir "/parset.cfg", parset);
+ Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+
+ MyGeometry::render();
+ MyGeometry::write();
+
+ // {{{ Set up grid
+ using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
+ auto grid = constructGrid<Grid>();
+
+ auto const refinements = parset.get<size_t>("grid.refinements");
+ grid->globalRefine(refinements);
+ size_t const fineVertexCount = grid->size(grid->maxLevel(), dims);
+
+ using GridView = Grid::LeafGridView;
+ GridView const leafView = grid->leafView();
+
+ MyFaces<GridView> myFaces(leafView);
+
+ // Neumann boundary
+ BoundaryPatch<GridView> const neumannBoundary(leafView);
+
+ // Frictional Boundary
+ BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
+ Dune::BitSetVector<1> frictionalNodes(fineVertexCount);
+ frictionalBoundary.getVertices(frictionalNodes);
+
+ // Dirichlet Boundary
+ Dune::BitSetVector<dims> noNodes(fineVertexCount);
+ Dune::BitSetVector<dims> dirichletNodes(fineVertexCount);
+ for (size_t i = 0; i < fineVertexCount; ++i) {
+ if (myFaces.right.containsVertex(i))
+ dirichletNodes[i][0] = true;
+
+ if (myFaces.lower.containsVertex(i))
+ dirichletNodes[i][1] = true;
+ }
+
+ // Set up functions for time-dependent boundary conditions
+ using Function = Dune::VirtualFunction<double, double>;
+ using FunctionMap = SharedPointerMap<std::string, Function>;
+ FunctionMap functions;
+ {
+ initPython();
+ Python::import("boundaryconditions")
+ .get("Functions")
+ .toC<typename FunctionMap::Base>(functions);
+ }
+ auto const &velocityDirichletFunction =
+ functions.get("velocityDirichletCondition"),
+ &neumannFunction = functions.get("neumannCondition");
+
+ using MyAssembler = MyAssembler<GridView, dims>;
+ using Matrix = MyAssembler::Matrix;
+ using LocalMatrix = Matrix::block_type;
+ using Vector = MyAssembler::Vector;
+ using LocalVector = Vector::block_type;
+ using ScalarMatrix = MyAssembler::ScalarMatrix;
+ using ScalarVector = MyAssembler::ScalarVector;
+ using LocalScalarVector = ScalarVector::block_type;
+
+ MyAssembler myAssembler(leafView);
+
+ MyBody<dims> const body(parset);
+
+ Matrix A, C, M;
+ myAssembler.assembleElasticity(body.getYoungModulus(),
+ body.getPoissonRatio(), A);
+ myAssembler.assembleViscosity(body.getShearViscosityField(),
+ body.getBulkViscosityField(), C);
+ myAssembler.assembleMass(body.getDensityField(), M);
+ EnergyNorm<Matrix, Vector> const ANorm(A), MNorm(M);
+ // Q: Does it make sense to weigh them in this manner?
+ SumNorm<Vector> const AMNorm(1.0, ANorm, 1.0, MNorm);
+
+ ScalarMatrix frictionalBoundaryMass;
+ myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
+ frictionalBoundaryMass);
+ EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
+ frictionalBoundaryMass);
+
+ // Assemble forces
+ Vector gravityFunctional;
+ myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional);
+
+ // Problem formulation: right-hand side
+ auto const computeExternalForces = [&](double _relativeTime, Vector &_ell) {
+ myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
+ _relativeTime);
+ _ell += gravityFunctional;
+ };
+ Vector ell(fineVertexCount);
+ computeExternalForces(0.0, ell);
+
+ // {{{ Initial conditions
+ using LinearFactory = SolverFactory<
+ dims, BlockNonlinearTNNMGProblem<ConvexProblem<
+ ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
+ Grid>;
+ ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
+
+ // TODO: clean up once generic lambdas arrive
+ auto const solveLinearProblem = [&](
+ Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
+ Vector const &_rhs, Vector &_x, EnergyNorm<Matrix, Vector> const &_norm,
+ Dune::ParameterTree const &_localParset) {
+ LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
+ refinements, *grid, _dirichletNodes);
+
+ typename LinearFactory::ConvexProblem convexProblem(
+ 1.0, _matrix, zeroNonlinearity, _rhs, _x);
+ typename LinearFactory::BlockProblem problem(parset, convexProblem);
+
+ auto multigridStep = factory.getSolver();
+ multigridStep->setProblem(_x, problem);
+ LoopSolver<Vector> solver(
+ multigridStep, _localParset.get<size_t>("maximumIterations"),
+ _localParset.get<double>("tolerance"), &_norm,
+ _localParset.get<Solver::VerbosityMode>("verbosity"),
+ false); // absolute error
+
+ solver.preprocess();
+ solver.solve();
+ };
+
+ // Solve the stationary problem
+ Vector u_initial(fineVertexCount);
+ u_initial = 0.0;
+ solveLinearProblem(dirichletNodes, A, ell, u_initial, ANorm,
+ parset.sub("u0.solver"));
+
+ ScalarVector alpha_initial(fineVertexCount);
+ alpha_initial = parset.get<double>("boundary.friction.initialLogState");
+ ScalarVector normalStress(fineVertexCount);
+ myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
+ body.getYoungModulus(),
+ body.getPoissonRatio(), u_initial);
+ MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"));
+ auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
+ frictionalBoundary, frictionInfo, normalStress);
+ myGlobalNonlinearity->updateLogState(alpha_initial);
+
+ Vector v_initial(fineVertexCount);
+ v_initial = 0.0;
+ {
+ double v_initial_const;
+ velocityDirichletFunction.evaluate(0.0, v_initial_const);
+ assert(v_initial_const == 0.0);
+ }
+
+ Vector a_initial(fineVertexCount);
+ a_initial = 0.0;
+ {
+ // We solve Ma = ell - [Au + Cv + Psi(v)]
+ Vector accelerationRHS(fineVertexCount);
+ {
+ accelerationRHS = 0.0;
+ Arithmetic::addProduct(accelerationRHS, A, u_initial);
+ Arithmetic::addProduct(accelerationRHS, C, v_initial);
+ // NOTE: We assume differentiability of Psi at 0 here!
+ myGlobalNonlinearity->addGradient(v_initial, accelerationRHS);
+ accelerationRHS *= -1.0;
+ accelerationRHS += ell;
+ }
+ solveLinearProblem(noNodes, M, accelerationRHS, a_initial, MNorm,
+ parset.sub("a0.solver"));
+ }
+ // }}}
+
+ Vector vertexCoordinates(fineVertexCount);
+ {
+ Dune::MultipleCodimMultipleGeomTypeMapper<
+ GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
+ for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it) {
+ auto const geometry = it->geometry();
+ assert(geometry.corners() == 1);
+ vertexCoordinates[vertexMapper.map(*it)] = geometry.corner(0);
+ }
+ }
+ FrictionWriter<ScalarVector, Vector> frictionWriter(
+ vertexCoordinates, frictionalNodes,
+ [](LocalVector const &x) { return x[0]; });
+ auto const reportFriction = [&](Vector const &_u, Vector const &_v,
+ ScalarVector const &_alpha) {
+ ScalarVector c;
+ myGlobalNonlinearity->coefficientOfFriction(_v, c);
+ frictionWriter.writeKinetics(_u, _v);
+ frictionWriter.writeOther(c, _alpha);
+ };
+ reportFriction(u_initial, v_initial, alpha_initial);
+
+ MyVTKWriter<typename MyAssembler::VertexBasis,
+ typename MyAssembler::CellBasis> const
+ vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
+
+ if (parset.get<bool>("io.writeVTK")) {
+ ScalarVector stress;
+ myAssembler.assembleVonMisesStress(
+ body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress);
+ vtkWriter.write(0, u_initial, v_initial, alpha_initial, stress);
+ }
+
+ SpecialWriter<GridView> specialVelocityWriter("specialVelocities",
+ leafView);
+ SpecialWriter<GridView> specialDisplacementWriter("specialDisplacements",
+ leafView);
+
+ // Set up TNNMG solver
+ using NonlinearFactory =
+ SolverFactory<dims, MyBlockProblem<ConvexProblem<
+ GlobalNonlinearity<Matrix, Vector>, Matrix>>,
+ Grid>;
+ NonlinearFactory factory(parset.sub("solver.tnnmg"), refinements, *grid,
+ dirichletNodes);
+ auto multigridStep = factory.getSolver();
+
+ std::fstream iterationWriter("iterations", std::fstream::out),
+ relaxationWriter("relaxation", std::fstream::out);
+
+ auto timeSteppingScheme =
+ initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
+ velocityDirichletFunction, dirichletNodes, M, A, C,
+ u_initial, v_initial, a_initial);
+ auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
+ parset.get<Config::stateModel>("boundary.friction.stateModel"),
+ alpha_initial, frictionalNodes,
+ parset.get<double>("boundary.friction.L"));
+
+ Vector v = v_initial;
+ Vector v_m(fineVertexCount);
+ ScalarVector alpha(fineVertexCount);
+
+ auto const timeSteps = parset.get<size_t>("timeSteps.number"),
+ maximumStateFPI = parset.get<size_t>("v.fpi.maximumIterations"),
+ maximumIterations =
+ parset.get<size_t>("v.solver.maximumIterations");
+ auto const tau = parset.get<double>("problem.finalTime") / timeSteps,
+ tolerance = parset.get<double>("v.solver.tolerance"),
+ fixedPointTolerance = parset.get<double>("v.fpi.tolerance"),
+ relaxation = parset.get<double>("v.fpi.relaxation"),
+ requiredReduction =
+ parset.get<double>("v.fpi.requiredReduction");
+ auto const printProgress = parset.get<bool>("io.printProgress");
+ auto const verbosity =
+ parset.get<Solver::VerbosityMode>("v.solver.verbosity");
+ for (size_t timeStep = 1; timeStep <= timeSteps; ++timeStep) {
+ if (printProgress)
+ std::cout << std::setw(7) << timeStep << " " << std::flush;
+
+ stateUpdater->nextTimeStep();
+ timeSteppingScheme->nextTimeStep();
+
+ auto const relativeTime = double(timeStep) / double(timeSteps);
+ computeExternalForces(relativeTime, ell);
+
+ Matrix velocityMatrix;
+ Vector velocityRHS(fineVertexCount);
+ Vector velocityIterate(fineVertexCount);
+
+ stateUpdater->setup(tau);
+ timeSteppingScheme->setup(ell, tau, relativeTime, velocityRHS,
+ velocityIterate, velocityMatrix);
+
+ LoopSolver<Vector> velocityProblemSolver(multigridStep, maximumIterations,
+ tolerance, &AMNorm, verbosity,
+ false); // absolute error
+
+ size_t iterationCounter;
+ auto solveVelocityProblem = [&](Vector &_velocityIterate,
+ ScalarVector const &_alpha) {
+ myGlobalNonlinearity->updateLogState(_alpha);
+
+ // NIT: Do we really need to pass u here?
+ typename NonlinearFactory::ConvexProblem convexProblem(
+ 1.0, velocityMatrix, *myGlobalNonlinearity, velocityRHS,
+ _velocityIterate);
+ typename NonlinearFactory::BlockProblem velocityProblem(parset,
+ convexProblem);
+ multigridStep->setProblem(_velocityIterate, velocityProblem);
+
+ velocityProblemSolver.preprocess();
+ velocityProblemSolver.solve();
+ iterationCounter = velocityProblemSolver.getResult().iterations;
+ };
+
+ Vector u;
+ Vector v_saved;
+ ScalarVector alpha_saved;
+ double lastStateCorrection;
+ for (size_t stateFPI = 1; stateFPI <= maximumStateFPI; ++stateFPI) {
+ timeSteppingScheme->extractOldVelocity(v_m);
+ v_m *= 0.5;
+ Arithmetic::addProduct(v_m, 0.5, v);
+
+ stateUpdater->solve(v_m);
+ stateUpdater->extractLogState(alpha);
+
+ if (stateFPI == 1)
+ relaxationWriter << "N ";
+ else {
+ double const stateCorrection =
+ stateEnergyNorm.diff(alpha, alpha_saved);
+ if (stateFPI <= 2 // lastStateCorrection is only set for stateFPI > 2
+ or stateCorrection < requiredReduction * lastStateCorrection)
+ relaxationWriter << "N ";
+ else {
+ alpha *= (1.0 - relaxation);
+ Arithmetic::addProduct(alpha, relaxation, alpha_saved);
+ relaxationWriter << "Y ";
+ }
+ lastStateCorrection = stateCorrection;
+ }
+
+ solveVelocityProblem(velocityIterate, alpha);
+ timeSteppingScheme->postProcess(velocityIterate);
+ timeSteppingScheme->extractDisplacement(u);
+ timeSteppingScheme->extractVelocity(v);
+
+ iterationWriter << iterationCounter << " ";
+ if (printProgress)
+ std::cout << '.' << std::flush;
+
+ if (stateFPI > 1) {
+ double const velocityCorrection = AMNorm.diff(v_saved, v);
+ if (velocityCorrection < fixedPointTolerance)
+ break;
+ }
+ if (stateFPI == maximumStateFPI)
+ DUNE_THROW(Dune::Exception, "FPI failed to converge");
+
+ alpha_saved = alpha;
+ v_saved = v;
+ }
+ if (printProgress)
+ std::cout << std::endl;
+
+ reportFriction(u, v, alpha);
+ iterationWriter << std::endl;
+ relaxationWriter << std::endl;
+
+ {
+ BasisGridFunction<typename MyAssembler::VertexBasis, Vector> velocity(
+ myAssembler.vertexBasis, v);
+ BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
+ displacement(myAssembler.vertexBasis, u);
+ specialVelocityWriter.write(velocity);
+ specialDisplacementWriter.write(displacement);
+ }
+
+ if (parset.get<bool>("io.writeVTK")) {
+ ScalarVector stress;
+ myAssembler.assembleVonMisesStress(body.getYoungModulus(),
+ body.getPoissonRatio(), u, stress);
+ vtkWriter.write(timeStep, u, v, alpha, stress);
+ }
+ }
+ iterationWriter.close();
+ relaxationWriter.close();
+
+ Python::stop();
+ }
+ catch (Dune::Exception &e) {
+ Dune::derr << "Dune reported error: " << e << std::endl;
+ }
+ catch (std::exception &e) {
+ std::cerr << "Standard exception: " << e.what() << std::endl;
+ }
+}