From 09bbd1b5fd796a22036df2939e8eec392238fb07 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Fri, 16 Feb 2018 16:21:42 +0100
Subject: [PATCH] .
---
dune-tectonic.pc.in | 2 +-
dune.module | 2 +-
dune/tectonic/globalratestatefriction.hh | 6 +-
dune/tectonic/myblockproblem.hh | 2 +-
src/CMakeLists.txt | 28 +
src/explicitgrid.hh | 7 -
src/gridselector.hh | 30 -
src/hdf5/frictionalboundary-writer.cc | 14 +-
src/hdf5/frictionalboundary-writer_tmpl.cc | 11 -
src/hdf5/patchinfo-writer.cc | 3 +-
src/hdf5/restart-io.cc | 22 +-
src/hdf5/surface-writer.cc | 6 +-
src/matrices.hh | 15 +-
src/multi-body-problem-2D.cfg | 21 +
src/multi-body-problem-3D.cfg | 24 +
src/multi-body-problem-data/bc.hh | 18 +
src/multi-body-problem-data/cuboidgeometry.cc | 190 +++++++
src/multi-body-problem-data/cuboidgeometry.hh | 77 +++
src/multi-body-problem-data/geometry.tex | 68 +++
src/multi-body-problem-data/midpoint.hh | 12 +
src/multi-body-problem-data/mybody.hh | 55 ++
.../myglobalfrictiondata.hh | 42 ++
src/multi-body-problem-data/mygrids.cc | 250 +++++++++
src/multi-body-problem-data/mygrids.hh | 91 +++
src/multi-body-problem-data/mygrids_tmpl.cc | 22 +
src/multi-body-problem-data/patchfunction.hh | 31 ++
.../segmented-function.hh | 34 ++
src/multi-body-problem-data/weakpatch.hh | 34 ++
src/multi-body-problem.cc | 523 ++++++++++++++++++
src/multi-body-problem.cfg | 71 +++
src/one-body-problem.cc | 31 +-
src/program_state.hh | 108 ++--
src/spatial-solving/fixedpointiterator.cc | 306 +++++++++-
src/spatial-solving/fixedpointiterator.hh | 53 --
src/spatial-solving/solverfactory.cc | 59 --
src/spatial-solving/solverfactory.hh | 48 --
src/spatial-solving/solverfactory_tmpl.cc | 22 -
src/time-stepping/rate.cc | 14 +-
src/time-stepping/rate.hh | 10 +-
src/time-stepping/rate/backward_euler.cc | 142 ++---
src/time-stepping/rate/backward_euler.hh | 15 +-
src/time-stepping/rate/newmark.cc | 116 ++--
src/time-stepping/rate/newmark.hh | 15 +-
src/time-stepping/rate/rateupdater.cc | 24 +-
src/time-stepping/rate/rateupdater.hh | 26 +-
src/time-stepping/rate_tmpl.cc | 11 +-
src/time-stepping/state.cc | 5 +-
src/time-stepping/state.hh | 4 +-
.../state/ageinglawstateupdater.cc | 33 +-
.../state/ageinglawstateupdater.hh | 20 +-
.../state/sliplawstateupdater.cc | 35 +-
.../state/sliplawstateupdater.hh | 28 +-
src/time-stepping/state/stateupdater.hh | 4 +-
src/time-stepping/state_tmpl.cc | 4 +-
54 files changed, 2307 insertions(+), 537 deletions(-)
create mode 100644 src/multi-body-problem-2D.cfg
create mode 100644 src/multi-body-problem-3D.cfg
create mode 100644 src/multi-body-problem-data/bc.hh
create mode 100644 src/multi-body-problem-data/cuboidgeometry.cc
create mode 100644 src/multi-body-problem-data/cuboidgeometry.hh
create mode 100644 src/multi-body-problem-data/geometry.tex
create mode 100644 src/multi-body-problem-data/midpoint.hh
create mode 100644 src/multi-body-problem-data/mybody.hh
create mode 100644 src/multi-body-problem-data/myglobalfrictiondata.hh
create mode 100644 src/multi-body-problem-data/mygrids.cc
create mode 100644 src/multi-body-problem-data/mygrids.hh
create mode 100644 src/multi-body-problem-data/mygrids_tmpl.cc
create mode 100644 src/multi-body-problem-data/patchfunction.hh
create mode 100644 src/multi-body-problem-data/segmented-function.hh
create mode 100644 src/multi-body-problem-data/weakpatch.hh
create mode 100644 src/multi-body-problem.cc
create mode 100644 src/multi-body-problem.cfg
diff --git a/dune-tectonic.pc.in b/dune-tectonic.pc.in
index d9c759b8..82d6d483 100644
--- a/dune-tectonic.pc.in
+++ b/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-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Requires: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
Libs: -L${libdir}
Cflags: -I${includedir}
diff --git a/dune.module b/dune.module
index e5c50280..8e630071 100644
--- a/dune.module
+++ b/dune.module
@@ -5,4 +5,4 @@
Module: dune-tectonic
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
+Depends: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
diff --git a/dune/tectonic/globalratestatefriction.hh b/dune/tectonic/globalratestatefriction.hh
index d6c61470..e233523d 100644
--- a/dune/tectonic/globalratestatefriction.hh
+++ b/dune/tectonic/globalratestatefriction.hh
@@ -32,7 +32,7 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
: restrictions(), localToGlobal(), zeroFriction() {
auto const gridView = frictionalBoundary.gridView();
Dune::MultipleCodimMultipleGeomTypeMapper<
- GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
+ GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
for (auto it = gridView.template begin<block_size>();
it != gridView.template end<block_size>(); ++it) {
auto const i = vertexMapper.index(*it);
@@ -44,8 +44,8 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
restrictions.emplace_back(weights[i], weightedNormalStress[i],
frictionInfo(geoToPoint(it->geometry())));
}
- assert(restrictions.size() == frictionalBoundary.numVertices());
- assert(localToGlobal.size() == frictionalBoundary.numVertices());
+ assert(restrictions.size() == (size_t) frictionalBoundary.numVertices());
+ assert(localToGlobal.size() == (size_t) frictionalBoundary.numVertices());
}
void updateAlpha(ScalarVector const &alpha) override {
diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index be3130ae..78cf36f4 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -57,7 +57,7 @@ class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
: Base(parset, problem),
maxEigenvalues_(problem.f.size()),
- localBisection(0.0, 1.0, 0.0, true, 0.0) {
+ localBisection(0.0, 1.0, 0.0, 0.0) {
for (size_t i = 0; i < problem.f.size(); ++i) {
LocalVectorType eigenvalues;
Dune::FMatrixHelp::eigenValues(problem.A[i][i], eigenvalues);
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 384e9885..3bb47317 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -19,6 +19,29 @@ set(SW_SOURCE_FILES
time-stepping/state.cc
vtk.cc
)
+
+set(MSW_SOURCE_FILES
+ assemblers.cc
+ enumparser.cc
+ hdf5/frictionalboundary-writer.cc
+ hdf5/iteration-writer.cc
+ hdf5/patchinfo-writer.cc
+ hdf5/restart-io.cc
+ hdf5/surface-writer.cc
+ hdf5/time-writer.cc
+ multi-body-problem-data/cuboidgeometry.cc
+ multi-body-problem-data/mygrids.cc
+ multi-body-problem.cc
+ spatial-solving/fixedpointiterator.cc
+ spatial-solving/solverfactory.cc
+ time-stepping/adaptivetimestepper.cc
+ time-stepping/coupledtimestepper.cc
+ time-stepping/rate.cc
+ time-stepping/rate/rateupdater.cc
+ time-stepping/state.cc
+ vtk.cc
+)
+
set(UGW_SOURCE_FILES
assemblers.cc # FIXME
one-body-problem-data/mygrid.cc
@@ -28,15 +51,20 @@ set(UGW_SOURCE_FILES
foreach(_dim 2 3)
set(_sw_target one-body-problem-${_dim}D)
+ set(_msw_target multi-body-problem-${_dim}D)
set(_ugw_target uniform-grid-writer-${_dim}D)
add_executable(${_sw_target} ${SW_SOURCE_FILES})
+ add_executable(${_msw_target} ${MSW_SOURCE_FILES})
add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
add_dune_ug_flags(${_sw_target})
+ add_dune_ug_flags(${_msw_target})
add_dune_ug_flags(${_ugw_target})
add_dune_hdf5_flags(${_sw_target})
+ add_dune_hdf5_flags(${_msw_target})
set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+ set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
endforeach()
diff --git a/src/explicitgrid.hh b/src/explicitgrid.hh
index fcbac18c..e69de29b 100644
--- a/src/explicitgrid.hh
+++ b/src/explicitgrid.hh
@@ -1,7 +0,0 @@
-#ifndef SRC_EXPLICITGRID_HH
-#define SRC_EXPLICITGRID_HH
-
-#include "gridselector.hh"
-
-using GridView = Grid::LeafGridView;
-#endif
diff --git a/src/gridselector.hh b/src/gridselector.hh
index 94fb9d85..e69de29b 100644
--- a/src/gridselector.hh
+++ b/src/gridselector.hh
@@ -1,30 +0,0 @@
-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#define WANT_ALUGRID 0
-#define WANT_UG 1
-
-#define WANT_GRID WANT_UG
-
-#if WANT_GRID == WANT_ALUGRID
-
-#if !HAVE_ALUGRID
-#error ALUGRID was requested but not found
-#endif
-#include <dune/grid/alugrid.hh>
-using Grid = Dune::ALUGrid<MY_DIM, MY_DIM, Dune::simplex, Dune::nonconforming>;
-
-#elif WANT_GRID == WANT_UG
-
-#if !HAVE_UG
-#error UG was requested but not found
-#endif
-#include <dune/grid/uggrid.hh>
-using Grid = Dune::UGGrid<MY_DIM>;
-
-#else
-
-#error requested a grid that does not exist
-
-#endif
diff --git a/src/hdf5/frictionalboundary-writer.cc b/src/hdf5/frictionalboundary-writer.cc
index e78b70f3..08bd86da 100644
--- a/src/hdf5/frictionalboundary-writer.cc
+++ b/src/hdf5/frictionalboundary-writer.cc
@@ -35,23 +35,27 @@ template <class ProgramState, class GridView>
template <class Friction>
void FrictionalBoundaryWriter<ProgramState, GridView>::write(
ProgramState const &programState, Friction &friction) {
+ // TODO
auto const frictionalBoundaryDisplacements =
- restrictToSurface(programState.u, frictionalBoundary_);
+ restrictToSurface(programState.u[0], frictionalBoundary_);
addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
frictionalBoundaryDisplacements);
auto const frictionalBoundaryVelocities =
- restrictToSurface(programState.v, frictionalBoundary_);
+ // TODO
+ restrictToSurface(programState.v[0], frictionalBoundary_);
addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
frictionalBoundaryVelocities);
auto const frictionalBoundaryStates =
- restrictToSurface(programState.alpha, frictionalBoundary_);
+ // TODO
+ restrictToSurface(programState.alpha[0], frictionalBoundary_);
addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
frictionalBoundaryStates);
- friction.updateAlpha(programState.alpha);
- auto const c = friction.coefficientOfFriction(programState.v);
+ // TODO
+ friction.updateAlpha(programState.alpha[0]);
+ auto const c = friction.coefficientOfFriction(programState.v[0]);
auto const frictionalBoundaryCoefficient =
restrictToSurface(c, frictionalBoundary_);
addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,
diff --git a/src/hdf5/frictionalboundary-writer_tmpl.cc b/src/hdf5/frictionalboundary-writer_tmpl.cc
index e39ffe0f..e69de29b 100644
--- a/src/hdf5/frictionalboundary-writer_tmpl.cc
+++ b/src/hdf5/frictionalboundary-writer_tmpl.cc
@@ -1,11 +0,0 @@
-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-using MyFriction = GlobalFriction<Matrix, Vector>;
-
-template class FrictionalBoundaryWriter<MyProgramState, GridView>;
-template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
- MyProgramState const &programState, MyFriction &friction);
diff --git a/src/hdf5/patchinfo-writer.cc b/src/hdf5/patchinfo-writer.cc
index 2df15c5f..fc723fc7 100644
--- a/src/hdf5/patchinfo-writer.cc
+++ b/src/hdf5/patchinfo-writer.cc
@@ -87,7 +87,8 @@ PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
template <class ProgramState, class VertexBasis, class GridView>
void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
ProgramState const &programState) {
- BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
+ // TODO
+ BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[0]);
auto const gridVelocity = gridEvaluator_.evaluate(velocity);
addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
}
diff --git a/src/hdf5/restart-io.cc b/src/hdf5/restart-io.cc
index d12ca34e..c5724e55 100644
--- a/src/hdf5/restart-io.cc
+++ b/src/hdf5/restart-io.cc
@@ -19,12 +19,13 @@ RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, size_t vertexCount)
template <class ProgramState>
void RestartIO<ProgramState>::write(ProgramState const &programState) {
- addEntry(displacementWriter_, programState.timeStep, programState.u);
- addEntry(velocityWriter_, programState.timeStep, programState.v);
- addEntry(accelerationWriter_, programState.timeStep, programState.a);
- addEntry(stateWriter_, programState.timeStep, programState.alpha);
+ // TODO
+ addEntry(displacementWriter_, programState.timeStep, programState.u[0]);
+ addEntry(velocityWriter_, programState.timeStep, programState.v[0]);
+ addEntry(accelerationWriter_, programState.timeStep, programState.a[0]);
+ addEntry(stateWriter_, programState.timeStep, programState.alpha[0]);
addEntry(weightedNormalStressWriter_, programState.timeStep,
- programState.weightedNormalStress);
+ programState.weightedNormalStress[0]);
addEntry(relativeTimeWriter_, programState.timeStep,
programState.relativeTime);
addEntry(relativeTimeIncrementWriter_, programState.timeStep,
@@ -35,12 +36,13 @@ template <class ProgramState>
void RestartIO<ProgramState>::read(size_t timeStep,
ProgramState &programState) {
programState.timeStep = timeStep;
- readEntry(displacementWriter_, timeStep, programState.u);
- readEntry(velocityWriter_, timeStep, programState.v);
- readEntry(accelerationWriter_, timeStep, programState.a);
- readEntry(stateWriter_, timeStep, programState.alpha);
+ // TODO
+ readEntry(displacementWriter_, timeStep, programState.u[0]);
+ readEntry(velocityWriter_, timeStep, programState.v[0]);
+ readEntry(accelerationWriter_, timeStep, programState.a[0]);
+ readEntry(stateWriter_, timeStep, programState.alpha[0]);
readEntry(weightedNormalStressWriter_, timeStep,
- programState.weightedNormalStress);
+ programState.weightedNormalStress[0]);
readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
}
diff --git a/src/hdf5/surface-writer.cc b/src/hdf5/surface-writer.cc
index bdab1213..207522bd 100644
--- a/src/hdf5/surface-writer.cc
+++ b/src/hdf5/surface-writer.cc
@@ -24,11 +24,13 @@ SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
template <class ProgramState, class GridView>
void SurfaceWriter<ProgramState, GridView>::write(
ProgramState const &programState) {
- auto const surfaceDisplacements = restrictToSurface(programState.u, surface_);
+ // TODO
+ auto const surfaceDisplacements = restrictToSurface(programState.u[0], surface_);
addEntry(surfaceDisplacementWriter_, programState.timeStep,
surfaceDisplacements);
- auto const surfaceVelocities = restrictToSurface(programState.v, surface_);
+ // TODO
+ auto const surfaceVelocities = restrictToSurface(programState.v[0], surface_);
addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
}
diff --git a/src/matrices.hh b/src/matrices.hh
index fc51d472..4fcaabfb 100644
--- a/src/matrices.hh
+++ b/src/matrices.hh
@@ -1,9 +1,16 @@
#ifndef SRC_MATRICES_HH
#define SRC_MATRICES_HH
-template <class Matrix> struct Matrices {
- Matrix elasticity;
- Matrix damping;
- Matrix mass;
+template <class Matrix> class Matrices {
+public:
+ std::vector<std::shared_ptr<Matrix>> elasticity;
+ std::vector<std::shared_ptr<Matrix>> damping;
+ std::vector<std::shared_ptr<Matrix>> mass;
+
+ Matrices(size_t n) {
+ elasticity.resize(n);
+ damping.resize(n);
+ mass.resize(n);
+ }
};
#endif
diff --git a/src/multi-body-problem-2D.cfg b/src/multi-body-problem-2D.cfg
new file mode 100644
index 00000000..a61cfba2
--- /dev/null
+++ b/src/multi-body-problem-2D.cfg
@@ -0,0 +1,21 @@
+# -*- mode:conf -*-
+[boundary.friction]
+smallestDiameter= 2e-3 # [m]
+
+[timeSteps]
+refinementTolerance = 1e-5
+
+[u0.solver]
+tolerance = 1e-8
+
+[a0.solver]
+tolerance = 1e-8
+
+[v.solver]
+tolerance = 1e-8
+
+[v.fpi]
+tolerance = 1e-5
+
+[solver.tnnmg.linear]
+tolerance = 1e-10
diff --git a/src/multi-body-problem-3D.cfg b/src/multi-body-problem-3D.cfg
new file mode 100644
index 00000000..3ff0794d
--- /dev/null
+++ b/src/multi-body-problem-3D.cfg
@@ -0,0 +1,24 @@
+# -*- mode:conf -*-
+[boundary.friction]
+smallestDiameter= 2e-2 # [m]
+
+[boundary.friction.weakening]
+patchType = Trapezoidal
+
+[timeSteps]
+refinementTolerance = 1e-5
+
+[u0.solver]
+tolerance = 1e-6
+
+[a0.solver]
+tolerance = 1e-6
+
+[v.solver]
+tolerance = 1e-6
+
+[v.fpi]
+tolerance = 1e-5
+
+[solver.tnnmg.linear]
+tolerance = 1e-10
diff --git a/src/multi-body-problem-data/bc.hh b/src/multi-body-problem-data/bc.hh
new file mode 100644
index 00000000..7c29cf08
--- /dev/null
+++ b/src/multi-body-problem-data/bc.hh
@@ -0,0 +1,18 @@
+#ifndef SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_BC_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;
+ y = (relativeTime <= 0.1)
+ ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
+ : finalVelocity;
+ }
+};
+
+class NeumannCondition : public Dune::VirtualFunction<double, double> {
+ void evaluate(double const &relativeTime, double &y) const { y = 0.0; }
+};
+#endif
diff --git a/src/multi-body-problem-data/cuboidgeometry.cc b/src/multi-body-problem-data/cuboidgeometry.cc
new file mode 100644
index 00000000..aaf4dfe0
--- /dev/null
+++ b/src/multi-body-problem-data/cuboidgeometry.cc
@@ -0,0 +1,190 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <fstream>
+
+#ifdef HAVE_CAIROMM
+#include <cairomm/context.h>
+#include <cairomm/fontface.h>
+#include <cairomm/surface.h>
+#endif
+
+#include "cuboidgeometry.hh"
+
+/*
+const CuboidGeometry::LocalVector CuboidGeometry::lift(const double v0, const double v1) {
+ vc = 0;
+ vc[0] = vc2D[0];
+ vc[1] = vc2D[1];
+}
+
+const CuboidGeometry::LocalVector& CuboidGeometry::A() const { return A_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::B() const { return B_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::C() const { return C_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::D() const { return D_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::X() const { return X_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::Y() const { return Y_; }
+double CuboidGeometry::depth() const { return depth_; }
+
+void CuboidGeometry::setupWeak(const LocalVector& weakOrigin, const double weakLength) {
+ lift(weakOrigin, X_);
+ const LocalVector2D Y({X_[0]+weakLength, X_[1]});
+ lift(Y, Y_);
+}
+*/
+
+#if MY_DIM == 3
+CuboidGeometry::CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength, const double length, const double width, const double depth_):
+ length_(length*lengthScale),
+ width_(width*lengthScale),
+ A(origin),
+ B({origin[0]+length_, origin[1], 0}),
+ C({origin[0]+length_, origin[1]+width_, 0}),
+ D({origin[0], origin[1]+width_, 0}),
+ X(weakOrigin),
+ Y({X[0]+weakLength, X[1], 0}),
+ depth(depth_*lengthScale) {}
+#else
+CuboidGeometry::CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength, const double length, const double width):
+ length_(length*lengthScale),
+ width_(width*lengthScale),
+ A(origin),
+ B({origin[0]+length_, origin[1]}),
+ C({origin[0]+length_, origin[1]+width_}),
+ D({origin[0], origin[1]+width_}),
+ X(weakOrigin),
+ Y({X[0]+weakLength, X[1]}) {}
+#endif
+
+
+void CuboidGeometry::write() const {
+ std::fstream writer("geometry", std::fstream::out);
+ writer << "A = " << A << std::endl;
+ writer << "B = " << B << std::endl;
+ writer << "C = " << C << std::endl;
+ writer << "D = " << D << std::endl;
+ writer << "X = " << X << std::endl;
+ writer << "Y = " << Y << std::endl;
+}
+
+void CuboidGeometry::render() const {
+#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 = [&](LocalVector2D const &v) { cr->move_to(v[0], -v[1]); };
+ auto const lineTo = [&](LocalVector2D const &v) { cr->line_to(v[0], -v[1]); };
+
+ cr->scale(widthScale, heightScale);
+ cr->translate(0.1, 0.1);
+ cr->set_line_width(0.0025);
+
+ // triangle
+ {
+ moveTo(reference::A);
+ lineTo(reference::B);
+ lineTo(reference::C);
+ cr->close_path();
+ cr->stroke();
+ }
+
+ // dashed lines
+ {
+ cr->save();
+ std::vector<double> dashPattern = { 0.005 };
+ cr->set_dash(dashPattern, 0);
+ moveTo(reference::Z);
+ lineTo(reference::Y);
+ moveTo(reference::U);
+ lineTo(reference::X);
+ cr->stroke();
+ cr->restore();
+ }
+
+ // fill viscoelastic region
+ {
+ cr->save();
+ setRGBColor(0x0097E0);
+ moveTo(reference::B);
+ lineTo(reference::K);
+ lineTo(reference::M);
+ cr->fill();
+ cr->restore();
+ }
+
+ // mark weakening region
+ {
+ cr->save();
+ setRGBColor(0x7AD3FF);
+ cr->set_line_width(0.005);
+ moveTo(reference::X);
+ lineTo(reference::Y);
+ cr->stroke();
+ cr->restore();
+ }
+
+ // mark points
+ {
+ auto const drawCircle = [&](LocalVector2D 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(reference::A);
+ drawCircle(reference::B);
+ drawCircle(reference::C);
+ drawCircle(reference::Y);
+ drawCircle(reference::X);
+ drawCircle(reference::Z);
+ drawCircle(reference::U);
+ drawCircle(reference::K);
+ drawCircle(reference::M);
+ drawCircle(reference::G);
+ drawCircle(reference::H);
+ drawCircle(reference::J);
+ drawCircle(reference::I);
+ cr->restore();
+ }
+
+ // labels
+ {
+ auto const label = [&](LocalVector2D const &v, std::string l) {
+ 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(D, "D");
+ label(X, "X");
+ label(Y, "Y");
+
+ cr->restore();
+ }
+
+ surface->write_to_png(filename);
+#endif
+}
diff --git a/src/multi-body-problem-data/cuboidgeometry.hh b/src/multi-body-problem-data/cuboidgeometry.hh
new file mode 100644
index 00000000..24245448
--- /dev/null
+++ b/src/multi-body-problem-data/cuboidgeometry.hh
@@ -0,0 +1,77 @@
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
+
+#include <dune/common/fvector.hh>
+
+#include "midpoint.hh"
+
+class CuboidGeometry {
+ public:
+ typedef Dune::FieldVector<double, MY_DIM> LocalVector;
+
+ constexpr static double const lengthScale = 1.0; // scaling factor
+
+ private:
+ const double length_;
+ const double width_;
+ //const double weakLength_; // length of the weak zone
+
+ public:
+ const LocalVector A;
+ const LocalVector B;
+ const LocalVector C;
+ const LocalVector D;
+ const LocalVector X;
+ const LocalVector Y;
+
+#if MY_DIM == 3
+ const double depth;
+
+ CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength = 0.20, const double length = 1.00, const double width = 0.27, const double depth = 0.60);
+#else
+ CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength = 0.20, const double length = 1.00, const double width = 0.27);
+#endif
+
+ void write() const;
+
+ void render() const;
+};
+#endif
+
+ /* from Elias
+ double const s = 1.0; // scaling factor
+
+ double const rightLeg = 0.27 * s;
+ double const leftLeg = 1.00 * s;
+ double const leftAngle = atan(rightLeg / leftLeg);
+ double const viscoHeight = 0.06 * s; // Height of the viscous bottom layer
+ double const weakLen = 0.20 * s; // Length of the weak zone
+
+ double const zDistance = 0.35;
+
+ LocalVector2D const A = {0, 0};
+ LocalVector2D const B = {leftLeg, -rightLeg};
+ LocalVector2D const C = {leftLeg, 0};
+
+ LocalVector2D const Z = {zDistance * s, 0};
+ LocalVector2D const Y = {zDistance * s, -zDistance *s / leftLeg *rightLeg};
+ LocalVector2D const X = {Y[0] - weakLen * std::cos(leftAngle),
+ Y[1] + weakLen *std::sin(leftAngle)};
+
+ LocalVector2D const U = {X[0], 0};
+
+ LocalVector2D const K = {B[0] - leftLeg * viscoHeight / rightLeg,
+ B[1] + viscoHeight};
+ LocalVector2D const M = {B[0], B[1] + viscoHeight};
+
+ LocalVector2D const G = midPoint(A, X);
+ LocalVector2D const H = midPoint(X, Y);
+ LocalVector2D const J = midPoint(Y, B);
+
+ LocalVector2D const I = {Y[0] + G[0], Y[1] + G[1]};
+
+ LocalVector2D const zenith = {0, 1};
+
+ LocalMatrix2D const rotation = {{std::cos(leftAngle), -std::sin(leftAngle)},
+ {std::sin(leftAngle), std::cos(leftAngle)}};
+ */
diff --git a/src/multi-body-problem-data/geometry.tex b/src/multi-body-problem-data/geometry.tex
new file mode 100644
index 00000000..32d63b6e
--- /dev/null
+++ b/src/multi-body-problem-data/geometry.tex
@@ -0,0 +1,68 @@
+\documentclass[tikz]{minimal}
+
+\usepackage{tikz}
+\usetikzlibrary{calc}
+\usetikzlibrary{decorations.pathreplacing}
+
+\usepackage{siunitx}
+
+\begin{document}
+\pgfmathsetmacro{\rightleg}{0.27}
+\pgfmathsetmacro{\leftleg}{1.00}
+\pgfmathsetmacro{\leftangle}{atan(\rightleg/\leftleg)}
+\begin{tikzpicture}[scale=12, rotate=\leftangle]
+ \pgfmathsetmacro{\mysin}{sin(\leftangle)}
+ \pgfmathsetmacro{\mycos}{cos(\leftangle)}
+ \pgfmathsetmacro{\viscoheight}{0.06}
+ \pgfmathsetmacro{\Zx}{0.35}
+ \pgfmathsetmacro{\weaklen}{0.20}
+
+ \coordinate (A) at (0,0);
+ \node at (A) [left] {A};
+ \coordinate (B) at (\leftleg,-\rightleg);
+ \node at (B) [right] {B};
+ \coordinate (C) at (\leftleg,0);
+ \node at (C) [right] {C};
+
+ \draw (A) -- (B) -- (C) -- node [above=.5cm, sloped] {$\overline{AC}=\SI{100}{cm}$} (A);
+
+ \coordinate (Z) at (\Zx,0);
+ \node at (Z) [above] {Z};
+ \coordinate (Y) at ($(\Zx,-\Zx/\leftleg * \rightleg)$);
+ \node at (Y) [below] {Y};
+ \coordinate (X) at ($(Y) + (-\weaklen*\mycos,\weaklen*\mysin)$);
+ \node at (X) [below] {X};
+ \path let \p1 = (X) in coordinate (U) at ($(\x1, 0)$);
+ \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] {$\overline{XY}=\SI{20}{cm}$} (Y);
+ \draw[dashed] (Y) -- (Z);
+ \draw[dashed] (U) -- (X);
+
+ \coordinate (K) at ($(B) + (-\leftleg * \viscoheight / \rightleg,\viscoheight)$);
+ \node at (K) [below] {K};
+ \coordinate (M) at ($(B) + (0, \viscoheight)$);
+ \node at (M) [right] {M};
+ \path (C) -- node [right=.5cm] {$\overline{CM} = \SI{21}{cm}$} (M);
+
+ \path[fill=blue] (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};
+
+ \node[align=left] at (0.5,-0.225) {
+ $Z$: coast line\\
+ $\overline{XY}$: velocity weakening zone\\
+ $BKM$: visco-elastic domain};
+\end{tikzpicture}
+
+\end{document}
diff --git a/src/multi-body-problem-data/midpoint.hh b/src/multi-body-problem-data/midpoint.hh
new file mode 100644
index 00000000..9b12a540
--- /dev/null
+++ b/src/multi-body-problem-data/midpoint.hh
@@ -0,0 +1,12 @@
+#ifndef SRC_MIDPOINT_HH
+#define SRC_MIDPOINT_HH
+
+#include <dune/matrix-vector/axpy.hh>
+
+template <class Vector> Vector midPoint(Vector const &x, Vector const &y) {
+ Vector ret(0);
+ Dune::MatrixVector::addProduct(ret, 0.5, x);
+ Dune::MatrixVector::addProduct(ret, 0.5, y);
+ return ret;
+}
+#endif
diff --git a/src/multi-body-problem-data/mybody.hh b/src/multi-body-problem-data/mybody.hh
new file mode 100644
index 00000000..08b5d4b9
--- /dev/null
+++ b/src/multi-body-problem-data/mybody.hh
@@ -0,0 +1,55 @@
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYBODY_HH
+#define SRC_MULTI_BODY_PROBLEM_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 "cuboidgeometry.hh"
+#include "segmented-function.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_;
+ SegmentedFunction const shearViscosityField_;
+ SegmentedFunction const bulkViscosityField_;
+ SegmentedFunction const densityField_;
+ Gravity<dimension> const gravityField_;
+};
+#endif
diff --git a/src/multi-body-problem-data/myglobalfrictiondata.hh b/src/multi-body-problem-data/myglobalfrictiondata.hh
new file mode 100644
index 00000000..d92e7223
--- /dev/null
+++ b/src/multi-body-problem-data/myglobalfrictiondata.hh
@@ -0,0 +1,42 @@
+#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
+
+#include <dune/common/function.hh>
+
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "patchfunction.hh"
+
+template <class LocalVector>
+class MyGlobalFrictionData : public GlobalFrictionData<LocalVector::dimension> {
+private:
+ using typename GlobalFrictionData<LocalVector::dimension>::VirtualFunction;
+
+public:
+ MyGlobalFrictionData(Dune::ParameterTree const &parset,
+ ConvexPolyhedron<LocalVector> const &segment)
+ : 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"), segment),
+ b_(parset.get<double>("strengthening.b"),
+ parset.get<double>("weakening.b"), segment),
+ mu0_(parset.get<double>("mu0")) {}
+
+ double const &C() const override { return C_; }
+ double const &L() const override { return L_; }
+ double const &V0() const override { return V0_; }
+ VirtualFunction const &a() const override { return a_; }
+ VirtualFunction const &b() const override { return b_; }
+ double 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/multi-body-problem-data/mygrids.cc b/src/multi-body-problem-data/mygrids.cc
new file mode 100644
index 00000000..c7e9cd0f
--- /dev/null
+++ b/src/multi-body-problem-data/mygrids.cc
@@ -0,0 +1,250 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+#include "mygrids.hh"
+#include "midpoint.hh"
+#include "../diameter.hh"
+
+#if MY_DIM == 3
+SimplexManager::SimplexManager(unsigned int shift) : shift_(shift) {}
+#endif
+
+// back-to-front, front-to-back, front-to-back
+void SimplexManager::addFromVerticesFBB(unsigned int U, unsigned int V,
+ unsigned int W) {
+#if MY_DIM == 3
+ unsigned int const U2 = U + shift_;
+ unsigned int const V2 = V + shift_;
+ unsigned int const W2 = W + shift_;
+
+ simplices_.push_back({ U, V, W, U2 });
+ simplices_.push_back({ V, V2, W2, U2 });
+ simplices_.push_back({ W, W2, U2, V });
+#else
+ simplices_.push_back({ U, V, W });
+#endif
+}
+
+// back-to-front, back-to-front, front-to-back
+void SimplexManager::addFromVerticesFFB(unsigned int U, unsigned int V,
+ unsigned int W) {
+#if MY_DIM == 3
+ unsigned int const U2 = U + shift_;
+ unsigned int const V2 = V + shift_;
+ unsigned int const W2 = W + shift_;
+
+ simplices_.push_back({ U, V, W, U2 });
+ simplices_.push_back({ V, V2, W, U2 });
+ simplices_.push_back({ V2, W, U2, W2 });
+#else
+ simplices_.push_back({ U, V, W });
+#endif
+}
+
+auto SimplexManager::getSimplices() -> SimplexList const & {
+ return simplices_;
+}
+
+// Fix: 3D case (still Elias code)
+template <class Grid> GridsConstructor<Grid>::GridsConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &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.A;
+ const auto& B = cuboidGeometry.B;
+ const auto& C = cuboidGeometry.C;
+ const auto& D = cuboidGeometry.D;
+
+ 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());
+ }
+}
+
+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, CuboidGeometry const &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, CuboidGeometry const &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.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
+ }
+
+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"
diff --git a/src/multi-body-problem-data/mygrids.hh b/src/multi-body-problem-data/mygrids.hh
new file mode 100644
index 00000000..aa62f7f8
--- /dev/null
+++ b/src/multi-body-problem-data/mygrids.hh
@@ -0,0 +1,91 @@
+#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>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.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, CuboidGeometry const &cuboidGeometry_);
+
+private:
+ CuboidGeometry const &cuboidGeometry;
+
+ bool isClose(double a, double b) {
+ return std::abs(a - b) < 1e-14 * cuboidGeometry.lengthScale;
+ };
+
+ bool isClose2(double a, double b) {
+ return std::abs(a - b) <
+ 1e-14 * cuboidGeometry.lengthScale * cuboidGeometry.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);
+};
+
+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 GridsConstructor {
+public:
+ GridsConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries_);
+
+ std::vector<std::shared_ptr<Grid>>& getGrids();
+
+ template <class GridView>
+ MyFaces<GridView> constructFaces(GridView const &gridView, CuboidGeometry const &cuboidGeometry);
+
+private:
+ std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries;
+ std::vector<Dune::GridFactory<Grid>> gridFactories;
+ std::vector<std::shared_ptr<Grid>> grids;
+};
+
+double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale);
+
+template <class Grid, class LocalVector>
+void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+ double smallestDiameter, double lengthScale);
+
+#endif
+
+
diff --git a/src/multi-body-problem-data/mygrids_tmpl.cc b/src/multi-body-problem-data/mygrids_tmpl.cc
new file mode 100644
index 00000000..99ae7cd6
--- /dev/null
+++ b/src/multi-body-problem-data/mygrids_tmpl.cc
@@ -0,0 +1,22 @@
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+#include "cuboidgeometry.hh"
+
+template class GridsConstructor<Grid>;
+
+template struct MyFaces<GridView>;
+template struct MyFaces<LevelGridView>;
+
+template MyFaces<GridView> GridsConstructor<Grid>::constructFaces(
+ GridView const &gridView, CuboidGeometry const &CuboidGeometry_);
+
+template MyFaces<LevelGridView> GridsConstructor<Grid>::constructFaces(
+ LevelGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
+
+template void refine<Grid, LocalVector>(
+ Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+ double smallestDiameter, double lengthScale);
diff --git a/src/multi-body-problem-data/patchfunction.hh b/src/multi-body-problem-data/patchfunction.hh
new file mode 100644
index 00000000..72cdc867
--- /dev/null
+++ b/src/multi-body-problem-data/patchfunction.hh
@@ -0,0 +1,31 @@
+#ifndef SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+class PatchFunction
+ : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
+ Dune::FieldVector<double, 1>> {
+private:
+ using Polyhedron = ConvexPolyhedron<Dune::FieldVector<double, MY_DIM>>;
+
+ double const v1_;
+ double const v2_;
+ Polyhedron const &segment_;
+
+public:
+ PatchFunction(double v1, double v2, Polyhedron const &segment)
+ : v1_(v1), v2_(v2), segment_(segment) {}
+
+ void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
+ Dune::FieldVector<double, 1> &y) const {
+ y = distance(x, segment_, 1e-6 * MyGeometry::lengthScale) <= 1e-5 ? v2_
+ : v1_;
+ }
+};
+
+#endif
diff --git a/src/multi-body-problem-data/segmented-function.hh b/src/multi-body-problem-data/segmented-function.hh
new file mode 100644
index 00000000..cab2dfef
--- /dev/null
+++ b/src/multi-body-problem-data/segmented-function.hh
@@ -0,0 +1,34 @@
+#ifndef SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include "cuboidgeometry.hh"
+
+class SegmentedFunction
+ : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
+ Dune::FieldVector<double, 1>> {
+private:
+ bool liesBelow(Dune::FieldVector<double, MY_DIM> const &x,
+ Dune::FieldVector<double, MY_DIM> const &y,
+ Dune::FieldVector<double, MY_DIM> const &z) const {
+ return x[1] + (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1];
+ };
+ bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const {
+ return liesBelow(MyGeometry::K, MyGeometry::M, z);
+ };
+
+ double const _v1;
+ double const _v2;
+
+public:
+ SegmentedFunction(double v1, double v2) : _v1(v1), _v2(v2) {}
+
+ void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
+ Dune::FieldVector<double, 1> &y) const {
+ y = insideRegion2(x) ? _v2 : _v1;
+ }
+};
+#endif
diff --git a/src/multi-body-problem-data/weakpatch.hh b/src/multi-body-problem-data/weakpatch.hh
new file mode 100644
index 00000000..e48f5725
--- /dev/null
+++ b/src/multi-body-problem-data/weakpatch.hh
@@ -0,0 +1,34 @@
+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_WEAKPATCH_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_WEAKPATCH_HH
+
+#include "cuboidgeometry.hh"
+
+template <class LocalVector>
+ConvexPolyhedron<LocalVector> getWeakPatch(Dune::ParameterTree const &parset, CuboidGeometry const &cuboidGeometry) {
+ ConvexPolyhedron<LocalVector> weakPatch;
+#if MY_DIM == 3
+ weakPatch.vertices.resize(4);
+ weakPatch.vertices[0] = weakPatch.vertices[2] = cuboidGeometry.X;
+ weakPatch.vertices[1] = weakPatch.vertices[3] = cuboidGeometry.Y;
+ for (size_t k = 0; k < 2; ++k) {
+ weakPatch.vertices[k][2] = -cuboidGeometry.depth / 2.0;
+ weakPatch.vertices[k + 2][2] = cuboidGeometry.depth / 2.0;
+ }
+ switch (parset.get<Config::PatchType>("patchType")) {
+ case Config::Rectangular:
+ break;
+ case Config::Trapezoidal:
+ weakPatch.vertices[1][0] += 0.05 * cuboidGeometry.lengthScale;
+ weakPatch.vertices[3][0] -= 0.05 * cuboidGeometry.lengthScale;
+ break;
+ default:
+ assert(false);
+ }
+#else
+ weakPatch.vertices.resize(2);
+ weakPatch.vertices[0] = cuboidGeometry.X;
+ weakPatch.vertices[1] = cuboidGeometry.Y;
+#endif
+ return weakPatch;
+};
+#endif
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
new file mode 100644
index 00000000..7533b579
--- /dev/null
+++ b/src/multi-body-problem.cc
@@ -0,0 +1,523 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#include <atomic>
+#include <cmath>
+#include <csignal>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <memory>
+
+#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/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/formatstring.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+
+
+/*
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/solvers/solver.hh>*/
+#include <dune/tnnmg/problem-classes/convexproblem.hh>
+
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/assemblers/nbodyassembler.hh>
+//#include <dune/contact/assemblers/dualmortarcoupling.hh>
+#include <dune/contact/projections/normalprojection.hh>
+
+#include <dune/tectonic/geocoordinate.hh>
+#include <dune/tectonic/myblockproblem.hh>
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/fufem/hdf5/file.hh>
+
+#include "assemblers.hh"
+#include "diameter.hh"
+#include "enumparser.hh"
+#include "enums.hh"
+#include "gridselector.hh"
+#include "hdf5-writer.hh"
+#include "hdf5/restart-io.hh"
+#include "matrices.hh"
+#include "program_state.hh"
+#include "multi-body-problem-data/bc.hh"
+#include "multi-body-problem-data/mybody.hh"
+#include "multi-body-problem-data/cuboidgeometry.hh"
+#include "multi-body-problem-data/myglobalfrictiondata.hh"
+#include "multi-body-problem-data/mygrids.hh"
+#include "multi-body-problem-data/weakpatch.hh"
+#include "spatial-solving/solverfactory.hh"
+#include "time-stepping/adaptivetimestepper.hh"
+#include "time-stepping/rate.hh"
+#include "time-stepping/state.hh"
+#include "time-stepping/updaters.hh"
+#include "vtk.hh"
+
+// for getcwd
+#include <unistd.h>
+
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+#define USE_OLD_TNNMG
+
+size_t const dims = MY_DIM;
+
+template <class GridView>
+void writeToVTK(const GridView& gridView, const std::string path, const std::string name) {
+ Dune::VTKWriter<GridView> vtkwriter(gridView);
+
+ //std::ofstream lStream( "garbage.txt" );
+ std::streambuf* lBufferOld = std::cout.rdbuf();
+ //std::cout.rdbuf(lStream.rdbuf());
+
+ vtkwriter.pwrite(name, path, path);
+
+ std::cout.rdbuf( lBufferOld );
+}
+
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+ Dune::ParameterTree parset;
+ Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem.cfg", parset);
+ Dune::ParameterTreeParser::readINITree(
+ Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem-%dD.cfg", dims), parset);
+ Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+ return parset;
+}
+
+static std::atomic<bool> terminationRequested(false);
+void handleSignal(int signum) { terminationRequested = true; }
+
+int main(int argc, char *argv[]) {
+ try {
+ Dune::MPIHelper::instance(argc, argv);
+
+ char buffer[256];
+ char *val = getcwd(buffer, sizeof(buffer));
+ if (val) {
+ std::cout << buffer << std::endl;
+ std::cout << argv[0] << std::endl;
+ }
+
+ auto const parset = getParameters(argc, argv);
+
+ using MyAssembler = MyAssembler<LeafGridView, dims>;
+ using Matrix = MyAssembler::Matrix;
+ using Vector = MyAssembler::Vector;
+ using LocalVector = Vector::block_type;
+ using ScalarMatrix = MyAssembler::ScalarMatrix;
+ using ScalarVector = MyAssembler::ScalarVector;
+ using field_type = Matrix::field_type;
+
+ // set up material properties of all bodies
+ MyBody<dims> const body(parset);
+
+ // set up cuboid geometries
+ const size_t bodyCount = parset.get<int>("problem.bodyCount");
+ std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries(bodyCount);
+
+ double const length = 1.00;
+ double const width = 0.27;
+ double const weakLength = 0.20;
+
+#if MY_DIM == 3
+ double const depth = 0.60;
+ // TODO: replace with make_shared
+ cuboidGeometries[0] = std::shared_ptr<CuboidGeometry>(new CuboidGeometry({0,0,0}, {0.2,width,0}, weakLength, length, width, depth));
+ for (size_t i=1; i<bodyCount; i++) {
+ cuboidGeometries[i] = std::shared_ptr<CuboidGeometry>(new CuboidGeometry(cuboidGeometries[i-1]->D, {0.6,i*width,0}, weakLength, length, width, depth));
+ }
+#else
+ // TODO: replace with make_shared
+ cuboidGeometries[0] = std::shared_ptr<CuboidGeometry>(new CuboidGeometry({0,0}, {0.2,width}, weakLength, length, width));
+ for (size_t i=1; i<bodyCount; i++) {
+ cuboidGeometries[i] = std::shared_ptr<CuboidGeometry>(new CuboidGeometry(cuboidGeometries[i-1]->D, {0.6,i*width}, weakLength, length, width));
+ }
+#endif
+
+ // set up (deformed) grids
+ GridsConstructor<Grid> gridConstructor(cuboidGeometries);
+ auto& grids = gridConstructor.getGrids();
+
+ using DeformedGridComplex = DeformedContinuaComplex<Grid, Vector>;
+ using DeformedGrid = DeformedGridComplex::DeformedGridType;
+ Dune::Contact::DeformedContinuaComplex<Grid, Vector> deformedGrids(grids);
+
+ using LeafGridView = DeformedGrid::LeafGridView;
+ using LevelGridView = DeformedGrid::LevelGridView;
+ std::vector<std::shared_ptr<LeafGridView>> leafViews(deformedGrids.size());
+ std::vector<std::shared_ptr<LevelGridView>> levelViews(deformedGrids.size());
+ std::vector<int> leafVertexCounts(deformedGrids.size());
+
+ using LeafFaces = MyFaces<LeafGridView>;
+ using LevelFaces = MyFaces<LevelGridView>;
+ std::vector<std::shared_ptr<LeafFaces>> leafFaces(deformedGrids.size());
+ std::vector<std::shared_ptr<LevelFaces>> levelFaces(deformedGrids.size());
+
+ std::vector<ConvexPolyhedron<LocalVector>> weakPatches(deformedGrids.size());
+
+ for (size_t i=0; i<deformedGrids.size(); i++) {
+ auto const &cuboidGeometry = *cuboidGeometries[i];
+
+ // define weak patch and refine grid
+ weakPatches[i] = getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"), cuboidGeometry);
+ refine(*grids[i], weakPatch, parset.get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
+
+ writeToVTK(grids[i]->leafGridView(), "", "grid" + std::to_string(i));
+
+ // determine minDiameter and maxDiameter
+ double minDiameter = std::numeric_limits<double>::infinity();
+ double maxDiameter = 0.0;
+ for (auto &&e : elements(grids[i]->leafGridView())) {
+ auto const geometry = e.geometry();
+ auto const diam = diameter(geometry);
+ minDiameter = std::min(minDiameter, diam);
+ maxDiameter = std::max(maxDiameter, diam);
+ }
+ std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
+ std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
+
+ leafViews[i] = std::make_shared<LeafGridView>(grids[i]->leafGridView());
+ levelViews[i] = std::make_shared<LevelGridView>(grids[i]->levelGridView(0));
+
+ leafVertexCounts[i] = leafViews[i]->size(dims);
+
+ leafFaces[i] = std::make_shared<LeafFaces>(*leafViews[i], cuboidGeometry);
+ levelFaces[i] = std::make_shared<LevelFaces>(*levelViews[i], cuboidGeometry);
+ }
+
+ // set up contact couplings
+ using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
+ using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
+
+ using CouplingPair = Dune::Contact::CouplingPair<DeformedGrid, DeformedGrid, field_type>;
+ std::vector<std::shared_ptr<CouplingPair>> couplings(bodyCount-1);
+
+ for (size_t i=0; i<couplings.size(); i++) {
+ auto nonmortarPatch = std::make_shared<LevelBoundaryPatch>(levelFaces[i]->upper);
+ auto mortarPatch = std::make_shared<LevelBoundaryPatch>(levelFaces[i+1]->lower);
+
+ auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
+ std::shared_ptr<CouplingPair::BackEndType> backend = nullptr;
+ couplings[i]->set(i, i+1, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
+ }
+
+ // set up dune-contact nBodyAssembler
+ std::vector<const DeformedGrid*> const_grids(bodyCount);
+ for (size_t i=0; i<const_grids.size(); i++) {
+ const_grids[i] = grids[i].get();
+ }
+
+ Dune::Contact::NBodyAssembler<DeformedGrid, Vector> nBodyAssembler(bodyCount, bodyCount-1);
+ nBodyAssembler.setGrids(const_grids);
+
+ for (size_t i=0; i<couplings.size(); i++) {
+ nBodyAssembler.setCoupling(*couplings[i], i);
+ }
+
+ nBodyAssembler.assembleTransferOperator();
+ nBodyAssembler.assembleObstacle();
+
+ // set up boundary conditions
+ std::vector<BoundaryPatch<LeafGridView>> neumannBoundaries(bodyCount);
+ std::vector<BoundaryPatch<LeafGridView>> frictionalBoundaries(bodyCount);
+ std::vector<BoundaryPatch<LeafGridView>> surfaces(bodyCount);
+
+ // Dirichlet boundary
+ std::vector<Dune::BitSetVector<dims>> noNodes(bodyCount);
+ std::vector<Dune::BitSetVector<dims>> dirichletNodes(bodyCount);
+
+ // set up functions for time-dependent boundary conditions
+ using Function = Dune::VirtualFunction<double, double>;
+ const Function& velocityDirichletFunction = VelocityDirichletCondition();
+ const Function& neumannFunction = NeumannCondition();
+
+ for (size_t i=0; i<grids.size(); i++) {
+ const auto& leafVertexCount = leafVertexCounts[i];
+
+ std::cout << "Grid" << i << " Number of DOFs: " << leafVertexCount << std::endl;
+
+ // Neumann boundary
+ auto& neumannBoundary = neumannBoundaries[i];
+ neumannBoundary = BoundaryPatch<LeafGridView>(*leafViews[i]);
+
+ // friction boundary
+ auto& frictionalBoundary = frictionalBoundaries[i];
+ if (i==0) {
+ frictionalBoundary = BoundaryPatch<LeafGridView>(leafFaces[i]->upper);
+ } else if (i==bodyCount-1) {
+ frictionalBoundary = BoundaryPatch<LeafGridView>(leafFaces[i]->lower);
+ } else {
+ frictionalBoundary = BoundaryPatch<LeafGridView>(leafFaces[i]->lower);
+ frictionalBoundary.addPatch(BoundaryPatch<LeafGridView>(leafFaces[i]->upper));
+ }
+ //surfaces[i] = BoundaryPatch<GridView>(myFaces.upper);
+
+ // TODO: Dirichlet Boundary
+ noNodes[i] = Dune::BitSetVector<dims>(leafVertexCount);
+ dirichletNodes[i] = Dune::BitSetVector<dims>(leafVertexCount);
+ auto& gridDirichletNodes = dirichletNodes[i];
+ for (int j=0; j<leafVertexCount; j++) {
+ if (leafFaces[i]->right.containsVertex(j))
+ gridDirichletNodes[j][0] = true;
+
+ if (leafFaces[i]->lower.containsVertex(j))
+ gridDirichletNodes[j][0] = true;
+
+ #if MY_DIM == 3
+ if (leafFaces[i]->front.containsVertex(j) || leafFaces[i]->back.containsVertex(j))
+ gridDirichletNodes[j][2] = true;
+ #endif
+ }
+ }
+
+ // set up individual assemblers for each body, assemble problem (matrices, forces, rhs)
+ std::vector<std::shared_ptr<MyAssembler>> assemblers(bodyCount);
+ Matrices<Matrix> matrices(bodyCount);
+
+ std::vector<Vector> gravityFunctionals(bodyCount);
+ std::vector<std::function<void(double, Vector&)>> externalForces(bodyCount);
+
+ for (size_t i=0; i<assemblers.size(); i++) {
+ auto& assembler = assemblers[i];
+ assembler = std::make_shared<MyAssembler>(*leafViews[i]);
+
+ assembler->assembleElasticity(body.getYoungModulus(), body.getPoissonRatio(), *matrices.elasticity[i]);
+ assembler->assembleViscosity(body.getShearViscosityField(), body.getBulkViscosityField(), *matrices.damping[i]);
+ assembler->assembleMass(body.getDensityField(), *matrices.mass[i]);
+
+ ScalarMatrix relativeFrictionalBoundaryMass;
+ assembler->assembleFrictionalBoundaryMass(frictionalBoundaries[i], relativeFrictionalBoundaryMass);
+ relativeFrictionalBoundaryMass /= frictionalBoundaries[i].area();
+ EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(relativeFrictionalBoundaryMass);
+
+ // assemble forces
+ assembler->assembleBodyForce(body.getGravityField(), gravityFunctionals[i]);
+
+ // Problem formulation: right-hand side
+ externalForces[i] =
+ [&](double _relativeTime, Vector &_ell) {
+ assembler->assembleNeumann(neumannBoundaries[i], _ell, neumannFunction,
+ _relativeTime);
+ _ell += gravityFunctionals[i];
+ };
+ }
+
+ /* Jonny 2bcontact
+ // make dirichlet bitfields containing dirichlet information for both grids
+ int size = rhs[0].size() + rhs[1].size();
+
+ Dune::BitSetVector<dims> totalDirichletNodes(size);
+
+ for (size_t i=0; i<dirichletNodes[0].size(); i++)
+ for (int j=0; j<dims; j++)
+ totalDirichletNodes[i][j] = dirichletNodes[0][i][j];
+
+ int offset = rhs[0].size();
+ for (size_t i=0; i<dirichletNodes[1].size(); i++)
+ for (int j=0; j<dims; j++)
+ totalDirichletNodes[offset + i][j] = dirichletNodes[1][i][j];
+
+ // assemble separate linear elasticity problems
+ std::array<MatrixType,2> stiffnessMatrix;
+ std::array<const MatrixType*, 2> submat;
+
+ for (size_t i=0; i<2; i++) {
+ OperatorAssembler<P1Basis,P1Basis> globalAssembler(*p1Basis[i],*p1Basis[i]);
+
+ double s = (i==0) ? E0 : E1;
+ StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> localAssembler(s, nu);
+
+ globalAssembler.assemble(localAssembler, stiffnessMatrix[i]);
+
+ submat[i] = &stiffnessMatrix[i];
+ }
+
+ MatrixType bilinearForm;
+ contactAssembler.assembleJacobian(submat, bilinearForm);
+ */
+
+ using MyProgramState = ProgramState<Vector, ScalarVector>;
+ MyProgramState programState(leafVertexCounts);
+ auto const firstRestart = parset.get<size_t>("io.restarts.first");
+ auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
+ auto const writeRestarts = parset.get<bool>("io.restarts.write");
+ auto const writeData = parset.get<bool>("io.data.write");
+ bool const handleRestarts = writeRestarts or firstRestart > 0;
+
+
+ auto dataFile =
+ writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+
+ auto restartFile = handleRestarts
+ ? std::make_unique<HDF5::File>(
+ "restarts.h5",
+ writeRestarts ? HDF5::Access::READWRITE
+ : HDF5::Access::READONLY)
+ : nullptr;
+
+
+ /* auto restartIO = handleRestarts
+ ? std::make_unique<RestartIO<MyProgramState>>(
+ *restartFile, leafVertexCounts)
+ : nullptr;*/
+/*
+ if (firstRestart > 0) // automatically adjusts the time and timestep
+ restartIO->read(firstRestart, programState);
+ else
+ programState.setupInitialConditions(parset, nBodyAssembler, externalForces,
+ matrices, assemblers, dirichletNodes,
+ noNodes, frictionalBoundaries, body);
+
+ */
+
+ // assemble friction
+ std::vector<std::shared_ptr<MyGlobalFrictionData<LocalVector>>> frictionInfo(weakPatches.size());
+ std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>> globalFriction(weakPatches.size());
+ for (size_t i=0; i<frictionInfo.size(); i++) {
+ frictionInfo[i] = std::make_shared<MyGlobalFrictionData<LocalVector>>(parset.sub("boundary.friction"), weakPatches[i]);
+ globalFriction[i] = assemblers[i]->assembleFrictionNonlinearity(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), frictionalBoundary, frictionInfo, programState.weightedNormalStress);
+ globalFriction[i]->updateAlpha(programState.alpha[i]);
+ }
+
+
+ /*
+ Vector vertexCoordinates(leafVertexCount);
+ {
+ Dune::MultipleCodimMultipleGeomTypeMapper<
+ GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView, Dune::mcmgVertexLayout());
+ for (auto &&v : vertices(leafView))
+ vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
+ }
+
+ using MyVertexBasis = typename MyAssembler::VertexBasis;
+ auto dataWriter =
+ writeData ? std::make_unique<
+ HDF5Writer<MyProgramState, MyVertexBasis, GridView>>(
+ *dataFile, vertexCoordinates, myAssembler.vertexBasis,
+ surface, frictionalBoundary, weakPatch)
+ : nullptr;
+ MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
+ myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
+
+
+ IterationRegister iterationCount;
+ auto const report = [&](bool initial = false) {
+ if (writeData) {
+ dataWriter->reportSolution(programState, *myGlobalFriction);
+ if (!initial)
+ dataWriter->reportIterations(programState, iterationCount);
+ dataFile->flush();
+ }
+
+ if (writeRestarts and !initial and
+ programState.timeStep % restartSpacing == 0) {
+ restartIO->write(programState);
+ restartFile->flush();
+ }
+
+ if (parset.get<bool>("io.printProgress"))
+ std::cout << "timeStep = " << std::setw(6) << programState.timeStep
+ << ", time = " << std::setw(12) << programState.relativeTime
+ << ", tau = " << std::setw(12) << programState.relativeTau
+ << std::endl;
+
+ if (parset.get<bool>("io.vtk.write")) {
+ ScalarVector stress;
+ myAssembler.assembleVonMisesStress(body.getYoungModulus(),
+ body.getPoissonRatio(),
+ programState.u, stress);
+ vtkWriter.write(programState.timeStep, programState.u, programState.v,
+ programState.alpha, stress);
+ }
+ };
+ report(true);
+ */
+
+
+ // Set up TNNMG solver
+ using NonlinearFactory = SolverFactory<
+ dims,
+ MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
+ Grid>;
+ NonlinearFactory factory(parset.sub("solver.tnnmg"), *grid, dirichletNodes);
+
+ using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
+ StateUpdater<ScalarVector, Vector>>;
+ MyUpdater current(
+ initRateUpdater(parset.get<Config::scheme>("timeSteps.scheme"),
+ velocityDirichletFunction, dirichletNodes, matrices,
+ programState.u, programState.v, programState.a),
+ initStateUpdater<ScalarVector, Vector>(
+ parset.get<Config::stateModel>("boundary.friction.stateModel"),
+ programState.alpha, *frictionalBoundary.getVertices(),
+ parset.get<double>("boundary.friction.L"),
+ parset.get<double>("boundary.friction.V0")));
+
+ auto const refinementTolerance =
+ parset.get<double>("timeSteps.refinementTolerance");
+ auto const mustRefine = [&](MyUpdater &coarseUpdater,
+ MyUpdater &fineUpdater) {
+ ScalarVector coarseAlpha;
+ coarseUpdater.state_->extractAlpha(coarseAlpha);
+
+ ScalarVector fineAlpha;
+ fineUpdater.state_->extractAlpha(fineAlpha);
+
+ return stateEnergyNorm.diff(fineAlpha, coarseAlpha) > refinementTolerance;
+ };
+
+ std::signal(SIGXCPU, handleSignal);
+ std::signal(SIGINT, handleSignal);
+ std::signal(SIGTERM, handleSignal);
+
+ AdaptiveTimeStepper<NonlinearFactory, MyUpdater,
+ EnergyNorm<ScalarMatrix, ScalarVector>>
+ adaptiveTimeStepper(factory, parset, myGlobalFriction, current,
+ programState.relativeTime, programState.relativeTau,
+ computeExternalForces, stateEnergyNorm, mustRefine);
+ while (!adaptiveTimeStepper.reachedEnd()) {
+ programState.timeStep++;
+
+ iterationCount = adaptiveTimeStepper.advance();
+
+ programState.relativeTime = adaptiveTimeStepper.relativeTime_;
+ programState.relativeTau = adaptiveTimeStepper.relativeTau_;
+ current.rate_->extractDisplacement(programState.u);
+ current.rate_->extractVelocity(programState.v);
+ current.rate_->extractAcceleration(programState.a);
+ current.state_->extractAlpha(programState.alpha);
+
+ report();
+
+ if (terminationRequested) {
+ std::cerr << "Terminating prematurely" << std::endl;
+ break;
+ }
+ }
+
+ } catch (Dune::Exception &e) {
+ Dune::derr << "Dune reported error: " << e << std::endl;
+ } catch (std::exception &e) {
+ std::cerr << "Standard exception: " << e.what() << std::endl;
+ }
+}
diff --git a/src/multi-body-problem.cfg b/src/multi-body-problem.cfg
new file mode 100644
index 00000000..b8358f55
--- /dev/null
+++ b/src/multi-body-problem.cfg
@@ -0,0 +1,71 @@
+# -*- mode:conf -*-
+gravity = 9.81 # [m/s^2]
+
+[io]
+data.write = true
+printProgress = false
+restarts.first = 0
+restarts.spacing= 20
+restarts.write = true
+vtk.write = false
+
+[problem]
+finalTime = 1000 # [s]
+bodyCount = 2
+
+[body]
+bulkModulus = 0.5e5 # [Pa]
+poissonRatio = 0.3 # [1]
+[body.elastic]
+density = 900 # [kg/m^3]
+shearViscosity = 1e3 # [Pas]
+bulkViscosity = 1e3 # [Pas]
+[body.viscoelastic]
+density = 1000 # [kg/m^3]
+shearViscosity = 1e4 # [Pas]
+bulkViscosity = 1e4 # [Pas]
+
+[boundary.friction]
+C = 10 # [Pa]
+mu0 = 0.7 # [ ]
+V0 = 5e-5 # [m/s]
+L = 2.25e-5 # [m]
+initialAlpha = 0 # [ ]
+stateModel = AgeingLaw
+frictionModel = Regularised
+[boundary.friction.weakening]
+a = 0.002 # [ ]
+b = 0.017 # [ ]
+[boundary.friction.strengthening]
+a = 0.020 # [ ]
+b = 0.005 # [ ]
+
+[timeSteps]
+scheme = newmark
+
+[u0.solver]
+maximumIterations = 100000
+verbosity = quiet
+
+[a0.solver]
+maximumIterations = 100000
+verbosity = quiet
+
+[v.solver]
+maximumIterations = 100000
+verbosity = quiet
+
+[v.fpi]
+maximumIterations = 10000
+lambda = 0.5
+
+[solver.tnnmg.linear]
+maxiumumIterations = 100000
+pre = 3
+cycle = 1 # 1 = V, 2 = W, etc.
+post = 3
+
+[solver.tnnmg.main]
+pre = 1
+multi = 5 # number of multigrid steps
+post = 0
diff --git a/src/one-body-problem.cc b/src/one-body-problem.cc
index d4d5ed56..dedc9162 100644
--- a/src/one-body-problem.cc
+++ b/src/one-body-problem.cc
@@ -31,9 +31,14 @@
#include <dune/fufem/formatstring.hh>
#include <dune/solvers/norms/energynorm.hh>
+
+
+/*
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/solvers/solver.hh>
#include <dune/tnnmg/problem-classes/convexproblem.hh>
+*/
+
#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/myblockproblem.hh>
@@ -62,13 +67,18 @@
#include "time-stepping/updaters.hh"
#include "vtk.hh"
+// for getcwd
+#include <unistd.h>
+
+#define USE_OLD_TNNMG
+
size_t const dims = MY_DIM;
Dune::ParameterTree getParameters(int argc, char *argv[]) {
Dune::ParameterTree parset;
- Dune::ParameterTreeParser::readINITree("one-body-problem.cfg", parset);
+ Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/one-body-problem.cfg", parset);
Dune::ParameterTreeParser::readINITree(
- Dune::Fufem::formatString("one-body-problem-%dD.cfg", dims), parset);
+ Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/one-body-problem-%dD.cfg", dims), parset);
Dune::ParameterTreeParser::readOptions(argc, argv, parset);
return parset;
}
@@ -79,6 +89,14 @@ void handleSignal(int signum) { terminationRequested = true; }
int main(int argc, char *argv[]) {
try {
Dune::MPIHelper::instance(argc, argv);
+
+ char buffer[256];
+ char *val = getcwd(buffer, sizeof(buffer));
+ if (val) {
+ std::cout << buffer << std::endl;
+ std::cout << argv[0] << std::endl;
+ }
+
auto const parset = getParameters(argc, argv);
MyGeometry::render();
@@ -140,6 +158,7 @@ int main(int argc, char *argv[]) {
#endif
}
+
// Set up functions for time-dependent boundary conditions
using Function = Dune::VirtualFunction<double, double>;
Function const &velocityDirichletFunction = VelocityDirichletCondition();
@@ -184,8 +203,10 @@ int main(int argc, char *argv[]) {
auto const writeData = parset.get<bool>("io.data.write");
bool const handleRestarts = writeRestarts or firstRestart > 0;
+ /*
auto dataFile =
writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+
auto restartFile = handleRestarts
? std::make_unique<HDF5::File>(
"restarts.h5",
@@ -214,7 +235,7 @@ int main(int argc, char *argv[]) {
Vector vertexCoordinates(leafVertexCount);
{
Dune::MultipleCodimMultipleGeomTypeMapper<
- GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
+ GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView, Dune::mcmgVertexLayout());
for (auto &&v : vertices(leafView))
vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
}
@@ -229,6 +250,7 @@ int main(int argc, char *argv[]) {
MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
+
IterationRegister iterationCount;
auto const report = [&](bool initial = false) {
if (writeData) {
@@ -261,6 +283,7 @@ int main(int argc, char *argv[]) {
};
report(true);
+
// Set up TNNMG solver
using NonlinearFactory = SolverFactory<
dims,
@@ -321,6 +344,8 @@ int main(int argc, char *argv[]) {
break;
}
}
+
+ */
} catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl;
} catch (std::exception &e) {
diff --git a/src/program_state.hh b/src/program_state.hh
index 1a63502f..73f9bc2e 100644
--- a/src/program_state.hh
+++ b/src/program_state.hh
@@ -3,10 +3,14 @@
#include <dune/common/parametertree.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/fufem/boundarypatch.hh>
#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
#include <dune/tectonic/body.hh>
#include "assemblers.hh"
@@ -18,32 +22,45 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
using Vector = VectorTEMPLATE;
using ScalarVector = ScalarVectorTEMPLATE;
- ProgramState(int leafVertexCount)
- : u(leafVertexCount),
- v(leafVertexCount),
- a(leafVertexCount),
- alpha(leafVertexCount),
- weightedNormalStress(leafVertexCount) {}
+ ProgramState(const std::vector<int>& leafVertexCounts)
+ : bodyCount(leafVertexCounts.size()) {
+ u.resize(bodyCount);
+ v.resize(bodyCount);
+ a.resize(bodyCount);
+ alpha.resize(bodyCount);
+ weightedNormalStress.resize(bodyCount);
+
+ for (size_t i=0; i<bodyCount; i++) {
+ size_t leafVertexCount = leafVertexCounts[i];
+ u[i].resize(leafVertexCount);
+ v[i].resize(leafVertexCount);
+ a[i].resize(leafVertexCount);
+ alpha[i].resize(leafVertexCount);
+ weightedNormalStress[i].resize(leafVertexCount);
+ }
+ }
// Set up initial conditions
template <class Matrix, class GridView>
void setupInitialConditions(
- Dune::ParameterTree const &parset,
- std::function<void(double, Vector &)> externalForces,
- Matrices<Matrix> const matrices,
- MyAssembler<GridView, Vector::block_type::dimension> const &myAssembler,
- Dune::BitSetVector<Vector::block_type::dimension> const &dirichletNodes,
- Dune::BitSetVector<Vector::block_type::dimension> const &noNodes,
- BoundaryPatch<GridView> const &frictionalBoundary,
- Body<Vector::block_type::dimension> const &body) {
+ const Dune::ParameterTree& parset,
+ const Dune::Contact::NBodyAssembler<typename GridView::Grid, Vector>& nBodyAssembler,
+ std::vector<std::function<void(double, Vector &)>> externalForces,
+ const Matrices<Matrix>& matrices,
+ const std::vector<std::shared_ptr<MyAssembler<GridView, Vector::block_type::dimension>>>& assemblers,
+ const std::vector<Dune::BitSetVector<Vector::block_type::dimension>>& dirichletNodes,
+ const std::vector<Dune::BitSetVector<Vector::block_type::dimension>>& noNodes,
+ const std::vector<BoundaryPatch<GridView>>& frictionalBoundaries,
+ const Body<Vector::block_type::dimension>& body) {
using LocalVector = typename Vector::block_type;
using LocalMatrix = typename Matrix::block_type;
auto constexpr dims = LocalVector::dimension;
+ /*
// Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
- Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
+ Dune::BitSetVector<dims> const &_dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrix,
Vector const &_rhs, Vector &_x,
Dune::ParameterTree const &_localParset) {
@@ -54,7 +71,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
- myAssembler.gridView.grid(), _dirichletNodes);
+ assemblers.gridView.grid(), _dirichletNodes);
typename LinearFactory::ConvexProblem convexProblem(
1.0, _matrix, zeroNonlinearity, _rhs, _x);
@@ -62,7 +79,11 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
auto multigridStep = factory.getStep();
multigridStep->setProblem(_x, problem);
+ //multigridStep->setProblem(_x);
EnergyNorm<Matrix, Vector> const norm(_matrix);
+
+
+
LoopSolver<Vector> solver(
multigridStep.get(), _localParset.get<size_t>("maximumIterations"),
_localParset.get<double>("tolerance"), &norm,
@@ -71,17 +92,28 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
solver.preprocess();
solver.solve();
+
+ Vector totalX = multigridStep->getSol();
+
+ // cleanup
+ delete(multigridStep);
+
+ nBodyAssembler.postprocess(totalX, x);
};
+ */
timeStep = 0;
relativeTime = 0.0;
relativeTau = 1e-6;
- Vector ell0(u.size());
- externalForces(relativeTime, ell0);
+ std::vector<Vector> ell0(bodyCount);
+ for (size_t i=0; i<bodyCount; i++) {
+ // Initial velocity
+ v[i] = 0.0;
- // Initial velocity
- v = 0.0;
+ ell0[i].resize(u[i].size());
+ externalForces[i](relativeTime, ell0[i]);
+ }
// Initial displacement: Start from a situation of minimal stress,
// which is automatically attained in the case [v = 0 = a].
@@ -91,29 +123,35 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// Initial acceleration: Computed in agreement with Ma = ell0 - Au
// (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
- Vector accelerationRHS = ell0;
- Arithmetic::subtractProduct(accelerationRHS, matrices.elasticity, u);
- solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a,
- parset.sub("a0.solver"));
+ std::vector<Vector> accelerationRHS = ell0;
+ for (size_t i=0; i<bodyCount; i++) {
+ // Initial state
+ alpha[i] = parset.get<double>("boundary.friction.initialAlpha");
+
+ // Initial normal stress
+ assemblers[i]->assembleWeightedNormalStress(
+ frictionalBoundaries[i], weightedNormalStress[i], body.getYoungModulus(),
+ body.getPoissonRatio(), u[i]);
- // Initial state
- alpha = parset.get<double>("boundary.friction.initialAlpha");
+ Dune::MatrixVector::subtractProduct(accelerationRHS[i], *matrices.elasticity[i], u[i]);
+ }
- // Initial normal stress
- myAssembler.assembleWeightedNormalStress(
- frictionalBoundary, weightedNormalStress, body.getYoungModulus(),
- body.getPoissonRatio(), u);
+ solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a,
+ parset.sub("a0.solver"));
}
public:
- Vector u;
- Vector v;
- Vector a;
- ScalarVector alpha;
- ScalarVector weightedNormalStress;
+ std::vector<Vector> u;
+ std::vector<Vector> v;
+ std::vector<Vector> a;
+ std::vector<ScalarVector> alpha;
+ std::vector<ScalarVector> weightedNormalStress;
double relativeTime;
double relativeTau;
size_t timeStep;
+
+private:
+ const size_t bodyCount;
};
#endif
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index a8c03459..02322caf 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -8,6 +8,19 @@
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/contact/assemblers/nbodyassembler.hh>
+#include <dune/contact/common/dualbasisadapter.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dune/functions/gridfunctions/gridfunction.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+
#include "../enums.hh"
#include "../enumparser.hh"
@@ -19,11 +32,13 @@ void FixedPointIterationCounter::operator+=(
multigridIterations += other.multigridIterations;
}
-template <class Factory, class Updaters, class ErrorNorm>
-FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
+template <class DeformedGrid, class Factory, class Updaters, class ErrorNorm>
+FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::FixedPointIterator(
+ const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
Factory &factory, Dune::ParameterTree const &parset,
std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
- : step_(factory.getStep()),
+ : nBodyAssembler_(nBodyAssembler),
+ step_(factory.getStep()),
parset_(parset),
globalFriction_(globalFriction),
fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
@@ -34,9 +49,9 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
errorNorm_(errorNorm) {}
-template <class Factory, class Updaters, class ErrorNorm>
+template <class DeformedGrid,class Factory, class Updaters, class ErrorNorm>
FixedPointIterationCounter
-FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
+FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::run(
Updaters updaters, Matrix const &velocityMatrix, Vector const &velocityRHS,
Vector &velocityIterate) {
EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
@@ -46,7 +61,7 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
size_t fixedPointIteration;
size_t multigridIterations = 0;
- ScalarVector alpha;
+ std::vector<ScalarVector> alpha;
updaters.state_->extractAlpha(alpha);
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
++fixedPointIteration) {
@@ -56,15 +71,22 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
velocityRHS, velocityIterate);
BlockProblem velocityProblem(parset_, convexProblem);
step_->setProblem(velocityIterate, velocityProblem);
+ //step_->setProblem(velocityIterate);
velocityProblemSolver.preprocess();
velocityProblemSolver.solve();
multigridIterations += velocityProblemSolver.getResult().iterations;
- Vector v_m;
+ std::vector<Vector> v_m;
updaters.rate_->extractOldVelocity(v_m);
- v_m *= 1.0 - lambda_;
- Arithmetic::addProduct(v_m, lambda_, velocityIterate);
+
+ for (size_t i=0; i<v_m.size(); i++) {
+ v_m[i] *= 1.0 - lambda_;
+ Arithmetic::addProduct(v_m[i], lambda_, velocityIterate[i]);
+ }
+
+ // compute relative velocities on contact boundaries
+ relativeVelocities(v_m);
// solve a state problem
updaters.state_->solve(v_m);
@@ -97,4 +119,270 @@ std::ostream &operator<<(std::ostream &stream,
<< ")";
}
+template <class DeformedGrid, class Factory, class Updaters, class ErrorNorm>
+void FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::relativeVelocities(std::vector<Vector>& v_m) const {
+ // adaptation of DualMortarCoupling::setup()
+
+ const size_t dim = DeformedGrid::dimension;
+ typedef typename DeformedGrid::LeafGridView GridView;
+
+ //cache of local bases
+ typedef Dune::PQkLocalFiniteElementCache<typename DeformedGrid::ctype, field_type, dim,1> FiniteElementCache1;
+ FiniteElementCache1 cache1;
+
+ // cache for the dual functions on the boundary
+ using DualCache = Dune::Contact::DualBasisAdapter<GridView, field_type>;
+ std::unique_ptr<DualCache> dualCache;
+ dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView, field_type> >();
+
+ // define FE grid functions
+ std::vector<BasisGridFunction<> > gridFunctions(nBodyAssembler_.nGrids());
+ for (size_t i=0; i<gridFunctions.size(); i++) {
+
+ }
+
+ const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
+ for (size_t i=0; i<contactCouplings.size(); i++) {
+ auto contactCoupling = contactCouplings[i];
+ auto glue = contactCoupling->getGlue();
+
+ // loop over all intersections
+ for (const auto& rIs : intersections(glue)) {
+ const auto& inside = rIs.inside();
+
+ if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+ continue;
+
+ const auto& outside = rIs.outside();
+
+ // types of the elements supporting the boundary segments in question
+ Dune::GeometryType nonmortarEType = inside.type();
+ Dune::GeometryType mortarEType = outside.type();
+
+ const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(nonmortarEType);
+ const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(mortarEType);
+
+ int noOfMortarVec = targetRefElement.size(dim);
+
+ Dune::GeometryType nmFaceType = domainRefElement.type(rIs.indexInInside(),1);
+ Dune::GeometryType mFaceType = targetRefElement.type(rIs.indexInOutside(),1);
+
+ // Select a quadrature rule
+ // 2 in 2d and for integration over triangles in 3d. If one (or both) of the two faces involved
+ // are quadrilaterals, then the quad order has to be risen to 3 (4).
+ int quadOrder = 2 + (!nmFaceType.isSimplex()) + (!mFaceType.isSimplex());
+ const auto& quadRule = Dune::QuadratureRules<ctype, dim-1>::rule(rIs.type(), quadOrder);
+
+ const auto& mortarFiniteElement = cache1.get(mortarEType);
+ dualCache->bind(inside, rIs.indexInInside());
+
+ std::vector<Dune::FieldVector<field_type,1> > mortarQuadValues, dualQuadValues;
+
+ const auto& rGeom = rIs.geometry();
+ const auto& rGeomOutside = rIs.geometryOutside();
+ const auto& rGeomInInside = rIs.geometryInInside();
+ const auto& rGeomInOutside = rIs.geometryInOutside();
+
+ int nNonmortarFaceNodes = domainRefElement.size(rIs.indexInInside(),1,dim);
+ std::vector<int> nonmortarFaceNodes;
+ for (int i=0; i<nNonmortarFaceNodes; i++) {
+ int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
+ nonmortarFaceNodes.push_back(faceIdxi);
+ }
+
+ for (const auto& quadPt : quadRule) {
+
+ // compute integration element of overlap
+ ctype integrationElement = rGeom.integrationElement(quadPt.position());
+
+ // quadrature point positions on the reference element
+ Dune::FieldVector<ctype,dim> nonmortarQuadPos = rGeomInInside.global(quadPt.position());
+ Dune::FieldVector<ctype,dim> mortarQuadPos = rGeomInOutside.global(quadPt.position());
+
+ // The current quadrature point in world coordinates
+ Dune::FieldVector<field_type,dim> nonmortarQpWorld = rGeom.global(quadPt.position());
+ Dune::FieldVector<field_type,dim> mortarQpWorld = rGeomOutside.global(quadPt.position());;
+
+ // the gap direction (normal * gapValue)
+ Dune::FieldVector<field_type,dim> gapVector = mortarQpWorld - nonmortarQpWorld;
+
+ //evaluate all shapefunctions at the quadrature point
+ //nonmortarFiniteElement.localBasis().evaluateFunction(nonmortarQuadPos,nonmortarQuadValues);
+ mortarFiniteElement.localBasis().evaluateFunction(mortarQuadPos,mortarQuadValues);
+ dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
+
+ // loop over all Lagrange multiplier shape functions
+ for (int j=0; j<nNonmortarFaceNodes; j++) {
+
+ int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
+ int rowIdx = globalToLocal[globalDomainIdx];
+
+ weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
+ * dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[globalDomainIdx]);
+
+ // loop over all mortar shape functions
+ for (int k=0; k<noOfMortarVec; k++) {
+
+ int colIdx = indexSet1.subIndex(outside, k, dim);
+ if (!mortarBoundary_.containsVertex(colIdx))
+ continue;
+
+ // Integrate over the product of two shape functions
+ field_type mortarEntry = integrationElement* quadPt.weight()* dualQuadValues[nonmortarFaceNodes[j]]* mortarQuadValues[k];
+
+ Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
+
+ }
+
+ }
+
+ }
+
+
+
+
+
+
+ }
+
+
+
+
+ ///////////////////////////////////
+ // reducing nonmortar boundary
+ /////////////////////////////////
+
+ // Get all fine grid boundary segments that are totally covered by the grid-glue segments
+ typedef std::pair<int,int> Pair;
+ std::map<Pair,ctype> coveredArea, fullArea;
+
+ // initialize with area of boundary faces
+ for (const auto& bIt : nonmortarBoundary_) {
+ const Pair p(indexSet0.index(bIt.inside()),bIt.indexInInside());
+ fullArea[p] = bIt.geometry().volume();
+ coveredArea[p] = 0;
+ }
+
+ // sum up the remote intersection areas to find out which are totally covered
+ for (const auto& rIs : intersections(glue))
+ coveredArea[Pair(indexSet0.index(rIs.inside()),rIs.indexInInside())] += rIs.geometry().volume();
+
+ // add all fine grid faces that are totally covered by the contact mapping
+ for (const auto& bIt : nonmortarBoundary_) {
+ const auto& inside = bIt.inside();
+ if(coveredArea[Pair(indexSet0.index(inside),bIt.indexInInside())]/
+ fullArea[Pair(indexSet0.index(inside),bIt.indexInInside())] >= coveredArea_)
+ boundaryWithMapping.addFace(inside, bIt.indexInInside());
+ }
+
+ //writeBoundary(boundaryWithMapping,debugPath_ + "relevantNonmortar");
+
+
+ /** \todo replace by all fine grid segments which are totally covered by the RemoteIntersections. */
+ //for (const auto& rIs : intersections(glue))
+ // boundaryWithMapping.addFace(rIs.inside(),rIs.indexInInside());
+
+ printf("contact mapping could be built for %d of %d boundary segments.\n",
+ boundaryWithMapping.numFaces(), nonmortarBoundary_.numFaces());
+
+ nonmortarBoundary_ = boundaryWithMapping;
+ mortarBoundary_.setup(gridView1);
+ for (const auto& rIs : intersections(glue))
+ if (nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+ mortarBoundary_.addFace(rIs.outside(),rIs.indexInOutside());
+
+
+ // Assemble the diagonal matrix coupling the nonmortar side with the lagrange multiplyers there
+ assembleNonmortarLagrangeMatrix();
+
+ // The weak obstacle vector
+ weakObstacle_.resize(nonmortarBoundary_.numVertices());
+ weakObstacle_ = 0;
+
+ // ///////////////////////////////////////////////////////////
+ // Get the occupation structure for the mortar matrix
+ // ///////////////////////////////////////////////////////////
+
+ /** \todo Also restrict mortar indices and don't use the whole grid level. */
+ Dune::MatrixIndexSet mortarIndices(nonmortarBoundary_.numVertices(), grid1_->size(dim));
+
+ // Create mapping from the global set of block dofs to the ones on the contact boundary
+ std::vector<int> globalToLocal;
+ nonmortarBoundary_.makeGlobalToLocal(globalToLocal);
+
+ // loop over all intersections
+ for (const auto& rIs : intersections(glue)) {
+
+ if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+ continue;
+
+ const auto& inside = rIs.inside();
+ const auto& outside = rIs.outside();
+
+ const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+ const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
+
+ int nDomainVertices = domainRefElement.size(dim);
+ int nTargetVertices = targetRefElement.size(dim);
+
+ for (int j=0; j<nDomainVertices; j++) {
+
+ int localDomainIdx = globalToLocal[indexSet0.subIndex(inside,j,dim)];
+
+ // if the vertex is not contained in the restricted contact boundary then dismiss it
+ if (localDomainIdx == -1)
+ continue;
+
+ for (int k=0; k<nTargetVertices; k++) {
+ int globalTargetIdx = indexSet1.subIndex(outside,k,dim);
+ if (!mortarBoundary_.containsVertex(globalTargetIdx))
+ continue;
+
+ mortarIndices.add(localDomainIdx, globalTargetIdx);
+ }
+ }
+ }
+
+ mortarIndices.exportIdx(mortarLagrangeMatrix_);
+
+ // Clear it
+ mortarLagrangeMatrix_ = 0;
+
+
+ //cache of local bases
+ FiniteElementCache1 cache1;
+
+ std::unique_ptr<DualCache> dualCache;
+ dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView0, field_type> >();
+
+ std::vector<Dune::FieldVector<ctype,dim> > avNormals;
+ avNormals = nonmortarBoundary_.getNormals();
+
+
+ }
+
+ // ///////////////////////////////////////
+ // Compute M = D^{-1} \hat{M}
+ // ///////////////////////////////////////
+
+ Dune::BCRSMatrix<MatrixBlock>& M = mortarLagrangeMatrix_;
+ Dune::BDMatrix<MatrixBlock>& D = nonmortarLagrangeMatrix_;
+
+ // First compute D^{-1}
+ D.invert();
+
+ // Then the matrix product D^{-1} \hat{M}
+ for (auto rowIt = M.begin(); rowIt != M.end(); ++rowIt) {
+ const auto rowIndex = rowIt.index();
+ for (auto& entry : *rowIt)
+ entry.leftmultiply(D[rowIndex][rowIndex]);
+ }
+
+ // weakObstacles in transformed basis = D^{-1}*weakObstacle_
+ for(size_t rowIdx=0; rowIdx<weakObstacle_.size(); rowIdx++)
+ weakObstacle_[rowIdx] *= D[rowIdx][rowIdx][0][0];
+
+ gridGlueBackend_->clear();
+}
+
#include "fixedpointiterator_tmpl.cc"
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index ba31c81e..e69de29b 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -1,53 +0,0 @@
-#ifndef SRC_SPATIAL_SOLVING_FIXEDPOINTITERATOR_HH
-#define SRC_SPATIAL_SOLVING_FIXEDPOINTITERATOR_HH
-
-#include <memory>
-
-#include <dune/common/parametertree.hh>
-
-#include <dune/solvers/norms/norm.hh>
-#include <dune/solvers/solvers/solver.hh>
-
-struct FixedPointIterationCounter {
- size_t iterations = 0;
- size_t multigridIterations = 0;
-
- void operator+=(FixedPointIterationCounter const &other);
-};
-
-std::ostream &operator<<(std::ostream &stream,
- FixedPointIterationCounter const &fpic);
-
-template <class Factory, class Updaters, class ErrorNorm>
-class FixedPointIterator {
- using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
- using Vector = typename Factory::Vector;
- using Matrix = typename Factory::Matrix;
- using ConvexProblem = typename Factory::ConvexProblem;
- using BlockProblem = typename Factory::BlockProblem;
- using Nonlinearity = typename ConvexProblem::NonlinearityType;
-
-public:
- FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
- std::shared_ptr<Nonlinearity> globalFriction,
- ErrorNorm const &errorNorm_);
-
- FixedPointIterationCounter run(Updaters updaters,
- Matrix const &velocityMatrix,
- Vector const &velocityRHS,
- Vector &velocityIterate);
-
-private:
- std::shared_ptr<typename Factory::Step> step_;
- Dune::ParameterTree const &parset_;
- std::shared_ptr<Nonlinearity> globalFriction_;
-
- size_t fixedPointMaxIterations_;
- double fixedPointTolerance_;
- double lambda_;
- size_t velocityMaxIterations_;
- double velocityTolerance_;
- Solver::VerbosityMode verbosity_;
- ErrorNorm const &errorNorm_;
-};
-#endif
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 3242724a..e69de29b 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -1,59 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
-#include <dune/fufem/assemblers/transferoperatorassembler.hh>
-#include <dune/solvers/solvers/solver.hh>
-
-#include "solverfactory.hh"
-
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
- Dune::ParameterTree const &parset, Grid const &grid,
- Dune::BitSetVector<dim> const &ignoreNodes)
- : baseEnergyNorm(linearBaseSolverStep),
- linearBaseSolver(&linearBaseSolverStep,
- parset.get<size_t>("linear.maxiumumIterations"),
- parset.get<double>("linear.tolerance"), &baseEnergyNorm,
- Solver::QUIET),
- transferOperators(grid.maxLevel()),
- multigridStep(
- std::make_shared<Step>(linearIterationStep, nonlinearSmoother)) {
- // linear iteration step
- linearIterationStep.setMGType(parset.get<int>("linear.cycle"),
- parset.get<int>("linear.pre"),
- parset.get<int>("linear.post"));
- linearIterationStep.basesolver_ = &linearBaseSolver;
- linearIterationStep.setSmoother(&linearPresmoother, &linearPostsmoother);
-
- // transfer operators
- for (auto &&x : transferOperators)
- x = new CompressedMultigridTransfer<Vector>;
- TransferOperatorAssembler<Grid>(grid)
- .assembleOperatorPointerHierarchy(transferOperators);
- linearIterationStep.setTransferOperators(transferOperators);
-
- // tnnmg iteration step
- multigridStep->setSmoothingSteps(parset.get<int>("main.pre"),
- parset.get<int>("main.multi"),
- parset.get<int>("main.post"));
- multigridStep->ignoreNodes_ = &ignoreNodes;
-}
-
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::~SolverFactory() {
- for (auto &&x : transferOperators)
- delete x;
-}
-
-template <size_t dim, class BlockProblem, class Grid>
-auto SolverFactory<dim, BlockProblem, Grid>::getStep()
- -> std::shared_ptr<Step> {
- return multigridStep;
-}
-
-#include "solverfactory_tmpl.cc"
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 8f05a572..e69de29b 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -1,48 +0,0 @@
-#ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
-#define SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/parametertree.hh>
-
-#include <dune/solvers/iterationsteps/multigridstep.hh>
-#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
-#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
-#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
-
-template <size_t dim, class BlockProblemTEMPLATE, class Grid>
-class SolverFactory {
-public:
- using BlockProblem = BlockProblemTEMPLATE;
- using ConvexProblem = typename BlockProblem::ConvexProblemType;
- using Vector = typename BlockProblem::VectorType;
- using Matrix = typename BlockProblem::MatrixType;
-
-private:
- using NonlinearSmoother = GenericNonlinearGS<BlockProblem>;
-
-public:
- using Step =
- TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
-
- SolverFactory(Dune::ParameterTree const &parset, Grid const &grid,
- Dune::BitSetVector<dim> const &ignoreNodes);
-
- ~SolverFactory();
-
- std::shared_ptr<Step> getStep();
-
-private:
- TruncatedBlockGSStep<Matrix, Vector> linearBaseSolverStep;
- EnergyNorm<Matrix, Vector> baseEnergyNorm;
- LoopSolver<Vector> linearBaseSolver;
- TruncatedBlockGSStep<Matrix, Vector> linearPresmoother;
- TruncatedBlockGSStep<Matrix, Vector> linearPostsmoother;
- MultigridStep<Matrix, Vector> linearIterationStep;
- std::vector<CompressedMultigridTransfer<Vector> *> transferOperators;
- NonlinearSmoother nonlinearSmoother;
- std::shared_ptr<Step> multigridStep;
-};
-#endif
diff --git a/src/spatial-solving/solverfactory_tmpl.cc b/src/spatial-solving/solverfactory_tmpl.cc
index 506d9042..e69de29b 100644
--- a/src/spatial-solving/solverfactory_tmpl.cc
+++ b/src/spatial-solving/solverfactory_tmpl.cc
@@ -1,22 +0,0 @@
-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
-
-#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
-
-template class SolverFactory<
- MY_DIM,
- MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
- Grid>;
-template class SolverFactory<
- MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
- ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
- Grid>;
diff --git a/src/time-stepping/rate.cc b/src/time-stepping/rate.cc
index 4d60bd9d..fe366fe7 100644
--- a/src/time-stepping/rate.cc
+++ b/src/time-stepping/rate.cc
@@ -9,20 +9,22 @@
template <class Vector, class Matrix, class Function, int dimension>
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
initRateUpdater(Config::scheme scheme,
- Function const &velocityDirichletFunction,
- Dune::BitSetVector<dimension> const &velocityDirichletNodes,
- Matrices<Matrix> const &matrices, Vector const &u_initial,
- Vector const &v_initial, Vector const &a_initial) {
+ const std::vector<Function>& velocityDirichletFunctions,
+ const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
+ const Matrices<Matrix>& matrices,
+ const std::vector<Vector>& u_initial,
+ const std::vector<Vector>& v_initial,
+ const std::vector<Vector>& a_initial) {
switch (scheme) {
case Config::Newmark:
return std::make_shared<Newmark<Vector, Matrix, Function, dimension>>(
matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
- velocityDirichletFunction);
+ velocityDirichletFunctions);
case Config::BackwardEuler:
return std::make_shared<
BackwardEuler<Vector, Matrix, Function, dimension>>(
matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
- velocityDirichletFunction);
+ velocityDirichletFunctions);
default:
assert(false);
}
diff --git a/src/time-stepping/rate.hh b/src/time-stepping/rate.hh
index 66b3962b..60268240 100644
--- a/src/time-stepping/rate.hh
+++ b/src/time-stepping/rate.hh
@@ -9,8 +9,10 @@
template <class Vector, class Matrix, class Function, int dimension>
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
initRateUpdater(Config::scheme scheme,
- Function const &velocityDirichletFunction,
- Dune::BitSetVector<dimension> const &velocityDirichletNodes,
- Matrices<Matrix> const &matrices, Vector const &u_initial,
- Vector const &v_initial, Vector const &a_initial);
+ const std::vector<Function>& velocityDirichletFunctions,
+ const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
+ const Matrices<Matrix>& matrices,
+ const std::vector<Vector>& u_initial,
+ const std::vector<Vector>& v_initial,
+ const std::vector<Vector>& a_initial);
#endif
diff --git a/src/time-stepping/rate/backward_euler.cc b/src/time-stepping/rate/backward_euler.cc
index 6664f02b..2694a970 100644
--- a/src/time-stepping/rate/backward_euler.cc
+++ b/src/time-stepping/rate/backward_euler.cc
@@ -1,88 +1,100 @@
-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/istl/matrixindexset.hh>
#include "backward_euler.hh"
template <class Vector, class Matrix, class Function, size_t dim>
BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
- Matrices<Matrix> const &_matrices, Vector const &_u_initial,
- Vector const &_v_initial, Vector const &_a_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- Function const &_dirichletFunction)
+ const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+ const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+ const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+ const std::vector<Function>& _dirichletFunctions)
: RateUpdater<Vector, Matrix, Function, dim>(
_matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
- _dirichletFunction) {}
+ _dirichletFunctions) {}
template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::setup(
- Vector const &ell, double _tau, double relativeTime, Vector &rhs,
- Vector &iterate, Matrix &AM) {
- this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
- this->tau = _tau;
-
- /* We start out with the formulation
-
- M a + C v + A u = ell
-
- Backward Euler means
-
- a1 = 1.0/tau ( v1 - v0 )
- u1 = tau v1 + u0
-
- in summary, we get at time t=1
-
- M [1.0/tau ( v1 - v0 )] + C v1
- + A [tau v1 + u0] = ell
-
- or
-
- 1.0/tau M v1 + C v1 + tau A v1
- = [1.0/tau M + C + tau A] v1
- = ell + 1.0/tau M v0 - A u0
- */
-
- // set up LHS (for fixed tau, we'd only really have to do this once)
- {
- Dune::MatrixIndexSet indices(this->matrices.elasticity.N(),
- this->matrices.elasticity.M());
- indices.import(this->matrices.elasticity);
- indices.import(this->matrices.mass);
- indices.import(this->matrices.damping);
- indices.exportIdx(AM);
+void BackwardEuler<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& ell,
+ double _tau,
+ double relativeTime,
+ std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+ std::vector<Matrix>& AM) {
+ for (size_t i=0; i<this->u.size(); i++) {
+ this->dirichletFunctions[i].evaluate(relativeTime, this->dirichletValues[i]);
+ this->tau = _tau;
+
+ /* We start out with the formulation
+
+ M a + C v + A u = ell
+
+ Backward Euler means
+
+ a1 = 1.0/tau ( v1 - v0 )
+ u1 = tau v1 + u0
+
+ in summary, we get at time t=1
+
+ M [1.0/tau ( v1 - v0 )] + C v1
+ + A [tau v1 + u0] = ell
+
+ or
+
+ 1.0/tau M v1 + C v1 + tau A v1
+ = [1.0/tau M + C + tau A] v1
+ = ell + 1.0/tau M v0 - A u0
+ */
+
+ // set up LHS (for fixed tau, we'd only really have to do this once)
+ Matrix& LHS = AM[i];
+ {
+ Dune::MatrixIndexSet indices(this->matrices.elasticity[i]->N(),
+ this->matrices.elasticity[i]->M());
+ indices.import(*this->matrices.elasticity[i]);
+ indices.import(*this->matrices.mass[i]);
+ indices.import(*this->matrices.damping[i]);
+ indices.exportIdx(LHS);
+ }
+ LHS = 0.0;
+ Dune::MatrixVector::addProduct(LHS, 1.0 / this->tau, *this->matrices.mass[i]);
+ Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
+ Dune::MatrixVector::addProduct(LHS, this->tau, *this->matrices.elasticity[i]);
+
+ // set up RHS
+ {
+ Vector& rhss = rhs[i];
+ rhss = ell[i];
+ Dune::MatrixVector::addProduct(rhss, 1.0 / this->tau, *this->matrices.mass[i],
+ this->v_o[i]);
+ Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+ }
+
+ iterate = this->v_o;
+
+ const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
+ for (size_t k = 0; k < dirichletNodess.size(); ++k)
+ for (size_t j = 0; j < dim; ++j)
+ if (this->dirichletNodes[k][j])
+ iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
}
- AM = 0.0;
- Arithmetic::addProduct(AM, 1.0 / this->tau, this->matrices.mass);
- Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
- Arithmetic::addProduct(AM, this->tau, this->matrices.elasticity);
-
- // set up RHS
- {
- rhs = ell;
- Arithmetic::addProduct(rhs, 1.0 / this->tau, this->matrices.mass,
- this->v_o);
- Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
- }
-
- iterate = this->v_o;
-
- for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
- for (size_t j = 0; j < dim; ++j)
- if (this->dirichletNodes[i][j])
- iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
}
template <class Vector, class Matrix, class Function, size_t dim>
void BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
- Vector const &iterate) {
+ const std::vector<Vector>& iterate) {
this->postProcessCalled = true;
this->v = iterate;
this->u = this->u_o;
- Arithmetic::addProduct(this->u, this->tau, this->v);
- this->a = this->v;
- this->a -= this->v_o;
- this->a /= this->tau;
+ for (size_t i=0; i<this->u.size(); i++) {
+ Dune::MatrixVector::addProduct(this->u[i], this->tau, this->v[i]);
+
+ Vector& ai = this->a[i];
+ ai = this->v[i];
+ ai -= this->v_o[i];
+ ai /= this->tau;
+ }
}
template <class Vector, class Matrix, class Function, size_t dim>
diff --git a/src/time-stepping/rate/backward_euler.hh b/src/time-stepping/rate/backward_euler.hh
index a9f3b3e7..245bb74b 100644
--- a/src/time-stepping/rate/backward_euler.hh
+++ b/src/time-stepping/rate/backward_euler.hh
@@ -4,14 +4,15 @@
template <class Vector, class Matrix, class Function, size_t dim>
class BackwardEuler : public RateUpdater<Vector, Matrix, Function, dim> {
public:
- BackwardEuler(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
- Vector const &_v_initial, Vector const &_a_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- Function const &_dirichletFunction);
+ BackwardEuler(const Matrices<Matrix> &_matrices, const std::vector<Vector> &_u_initial,
+ const std::vector<Vector> &_v_initial, const std::vector<Vector> &_a_initial,
+ const std::vector<Dune::BitSetVector<dim> > &_dirichletNodes,
+ const std::vector<Function> &_dirichletFunctions);
- void setup(Vector const &, double, double, Vector &, Vector &,
- Matrix &) override;
- void postProcess(Vector const &) override;
+ void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
+ std::vector<Matrix>&) override;
+
+ void postProcess(const std::vector<Vector>&) override;
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
const override;
diff --git a/src/time-stepping/rate/newmark.cc b/src/time-stepping/rate/newmark.cc
index 99076858..53bf32c4 100644
--- a/src/time-stepping/rate/newmark.cc
+++ b/src/time-stepping/rate/newmark.cc
@@ -1,27 +1,29 @@
-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/istl/matrixindexset.hh>
#include "newmark.hh"
template <class Vector, class Matrix, class Function, size_t dim>
Newmark<Vector, Matrix, Function, dim>::Newmark(
- Matrices<Matrix> const &_matrices, Vector const &_u_initial,
- Vector const &_v_initial, Vector const &_a_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- Function const &_dirichletFunction)
+ const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+ const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+ const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+ const std::vector<Function>& _dirichletFunctions)
: RateUpdater<Vector, Matrix, Function, dim>(
_matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
- _dirichletFunction) {}
+ _dirichletFunctions) {}
template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
+void Newmark<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& ell,
double _tau,
double relativeTime,
- Vector &rhs, Vector &iterate,
- Matrix &AM) {
- this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
- this->tau = _tau;
+ std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+ std::vector<Matrix>& AM) {
+ for (size_t i=0; i<this->u.size(); i++) {
+ this->dirichletFunctions[i].evaluate(relativeTime, this->dirichletValues[i]);
+ this->tau = _tau;
- /* We start out with the formulation
+ /* We start out with the formulation
M a + C v + A u = ell
@@ -41,58 +43,64 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
= [2/tau M + C + tau/2 A] v1
= ell + 2/tau M v0 + M a0
- tau/2 A v0 - A u0
- */
-
- // set up LHS (for fixed tau, we'd only really have to do this once)
- {
- Dune::MatrixIndexSet indices(this->matrices.elasticity.N(),
- this->matrices.elasticity.M());
- indices.import(this->matrices.elasticity);
- indices.import(this->matrices.mass);
- indices.import(this->matrices.damping);
- indices.exportIdx(AM);
+ */
+
+ // set up LHS (for fixed tau, we'd only really have to do this once)
+ Matrix& LHS = AM[i];
+ {
+ Dune::MatrixIndexSet indices(this->matrices.elasticity[i]->N(),
+ this->matrices.elasticity[i]->M());
+ indices.import(*this->matrices.elasticity[i]);
+ indices.import(*this->matrices.mass[i]);
+ indices.import(*this->matrices.damping[i]);
+ indices.exportIdx(LHS);
+ }
+ LHS = 0.0;
+ Dune::MatrixVector::addProduct(LHS, 2.0 / this->tau, *this->matrices.mass[i]);
+ Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
+ Dune::MatrixVector::addProduct(LHS, this->tau / 2.0, *this->matrices.elasticity[i]);
+
+ // set up RHS
+ {
+ Vector& rhss = rhs[i];
+ rhss = ell[i];
+ Dune::MatrixVector::addProduct(rhss, 2.0 / this->tau, *this->matrices.mass[i],
+ this->v_o[i]);
+ Dune::MatrixVector::addProduct(rhss, *this->matrices.mass[i], this->a_o[i]);
+ Dune::MatrixVector::subtractProduct(rhss, this->tau / 2.0, *this->matrices.elasticity[i],
+ this->v_o[i]);
+ Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+ }
+
+ iterate = this->v_o;
+
+ const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
+ for (size_t k = 0; k < dirichletNodess.size(); ++k)
+ for (size_t j = 0; j < dim; ++j)
+ if (this->dirichletNodess[k][j])
+ iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
}
- AM = 0.0;
- Arithmetic::addProduct(AM, 2.0 / this->tau, this->matrices.mass);
- Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
- Arithmetic::addProduct(AM, this->tau / 2.0, this->matrices.elasticity);
-
- // set up RHS
- {
- rhs = ell;
- Arithmetic::addProduct(rhs, 2.0 / this->tau, this->matrices.mass,
- this->v_o);
- Arithmetic::addProduct(rhs, this->matrices.mass, this->a_o);
- Arithmetic::subtractProduct(rhs, this->tau / 2.0, this->matrices.elasticity,
- this->v_o);
- Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
- }
-
- iterate = this->v_o;
-
- for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
- for (size_t j = 0; j < dim; ++j)
- if (this->dirichletNodes[i][j])
- iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
}
template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::postProcess(
- Vector const &iterate) {
+void Newmark<Vector, Matrix, Function, dim>::postProcess(const std::vector<Vector>& iterate) {
this->postProcessCalled = true;
this->v = iterate;
// u1 = tau/2 ( v1 + v0 ) + u0
this->u = this->u_o;
- Arithmetic::addProduct(this->u, this->tau / 2.0, this->v);
- Arithmetic::addProduct(this->u, this->tau / 2.0, this->v_o);
-
- // a1 = 2/tau ( v1 - v0 ) - a0
- this->a = 0.0;
- Arithmetic::addProduct(this->a, 2.0 / this->tau, this->v);
- Arithmetic::subtractProduct(this->a, 2.0 / this->tau, this->v_o);
- Arithmetic::subtractProduct(this->a, 1.0, this->a_o);
+
+ for (size_t i=0; i<this->u.size(); i++) {
+ Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v[i]);
+ Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v_o[i]);
+
+ // a1 = 2/tau ( v1 - v0 ) - a0
+ this->a[i] = 0.0;
+ Dune::MatrixVector::addProduct(this->a[i], 2.0 / this->tau, this->v[i]);
+ Dune::MatrixVector::subtractProduct(this->a[i], 2.0 / this->tau, this->v_o[i]);
+ Dune::MatrixVector::subtractProduct(this->a[i], 1.0, this->a_o[i]);
+ }
}
template <class Vector, class Matrix, class Function, size_t dim>
diff --git a/src/time-stepping/rate/newmark.hh b/src/time-stepping/rate/newmark.hh
index 400eff03..be5cb844 100644
--- a/src/time-stepping/rate/newmark.hh
+++ b/src/time-stepping/rate/newmark.hh
@@ -4,14 +4,15 @@
template <class Vector, class Matrix, class Function, size_t dim>
class Newmark : public RateUpdater<Vector, Matrix, Function, dim> {
public:
- Newmark(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
- Vector const &_v_initial, Vector const &_a_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- Function const &_dirichletFunction);
+ Newmark(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+ const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+ const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+ const std::vector<Function>& _dirichletFunctions);
- void setup(Vector const &, double, double, Vector &, Vector &,
- Matrix &) override;
- void postProcess(Vector const &) override;
+ void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
+ std::vector<Matrix>&) override;
+
+ void postProcess(const std::vector<Vector>&) override;
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
const override;
diff --git a/src/time-stepping/rate/rateupdater.cc b/src/time-stepping/rate/rateupdater.cc
index 765b6b64..23d09536 100644
--- a/src/time-stepping/rate/rateupdater.cc
+++ b/src/time-stepping/rate/rateupdater.cc
@@ -6,16 +6,16 @@
template <class Vector, class Matrix, class Function, size_t dim>
RateUpdater<Vector, Matrix, Function, dim>::RateUpdater(
- Matrices<Matrix> const &_matrices, Vector const &_u_initial,
- Vector const &_v_initial, Vector const &_a_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- Function const &_dirichletFunction)
+ const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+ const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+ const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+ const std::vector<Function>& _dirichletFunctions)
: matrices(_matrices),
u(_u_initial),
v(_v_initial),
a(_a_initial),
dirichletNodes(_dirichletNodes),
- dirichletFunction(_dirichletFunction) {}
+ dirichletFunctions(_dirichletFunctions) {}
template <class Vector, class Matrix, class Function, size_t dim>
void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() {
@@ -26,17 +26,15 @@ void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() {
}
template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractDisplacement(
- Vector &displacement) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractDisplacement(std::vector<Vector>& displacements) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
- displacement = u;
+ displacements = u;
}
template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(
- Vector &velocity) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(std::vector<Vector>& velocity) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
@@ -44,14 +42,12 @@ void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(
}
template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractOldVelocity(
- Vector &oldVelocity) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractOldVelocity(std::vector<Vector>& oldVelocity) const {
oldVelocity = v_o;
}
template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractAcceleration(
- Vector &acceleration) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractAcceleration(std::vector<Vector>& acceleration) const {
if (!postProcessCalled)
DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
diff --git a/src/time-stepping/rate/rateupdater.hh b/src/time-stepping/rate/rateupdater.hh
index 7ba9c7ba..6e2a9863 100644
--- a/src/time-stepping/rate/rateupdater.hh
+++ b/src/time-stepping/rate/rateupdater.hh
@@ -10,10 +10,10 @@
template <class Vector, class Matrix, class Function, size_t dim>
class RateUpdater {
protected:
- RateUpdater(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
- Vector const &_v_initial, Vector const &_a_initial,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- Function const &_dirichletFunction);
+ RateUpdater(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+ const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+ const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+ const std::vector<Function>& _dirichletFunctions);
public:
void nextTimeStep();
@@ -21,22 +21,22 @@ class RateUpdater {
Vector &rhs, Vector &iterate, Matrix &AB) = 0;
void virtual postProcess(Vector const &iterate) = 0;
- void extractDisplacement(Vector &displacement) const;
- void extractVelocity(Vector &velocity) const;
- void extractOldVelocity(Vector &velocity) const;
- void extractAcceleration(Vector &acceleration) const;
+ void extractDisplacement(std::vector<Vector>& displacements) const;
+ void extractVelocity(std::vector<Vector>& velocity) const;
+ void extractOldVelocity(std::vector<Vector>& velocity) const;
+ void extractAcceleration(std::vector<Vector>& acceleration) const;
std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> virtual clone()
const = 0;
protected:
- Matrices<Matrix> const &matrices;
- Vector u, v, a;
- Dune::BitSetVector<dim> const &dirichletNodes;
- Function const &dirichletFunction;
+ const Matrices<Matrix>& matrices;
+ std::vector<Vector> u, v, a;
+ const std::vector<Dune::BitSetVector<dim>>& dirichletNodes;
+ const std::vector<Function>& dirichletFunctions;
double dirichletValue;
- Vector u_o, v_o, a_o;
+ std::vector<Vector> u_o, v_o, a_o;
double tau;
bool postProcessCalled = true;
diff --git a/src/time-stepping/rate_tmpl.cc b/src/time-stepping/rate_tmpl.cc
index 728d55c4..92c15bdd 100644
--- a/src/time-stepping/rate_tmpl.cc
+++ b/src/time-stepping/rate_tmpl.cc
@@ -6,7 +6,10 @@ using Function = Dune::VirtualFunction<double, double>;
template std::shared_ptr<RateUpdater<Vector, Matrix, Function, MY_DIM>>
initRateUpdater<Vector, Matrix, Function, MY_DIM>(
- Config::scheme scheme, Function const &velocityDirichletFunction,
- Dune::BitSetVector<MY_DIM> const &velocityDirichletNodes,
- Matrices<Matrix> const &matrices, Vector const &u_initial,
- Vector const &v_initial, Vector const &a_initial);
+ Config::scheme scheme,
+ const std::vector<Function>& velocityDirichletFunctions,
+ const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
+ const Matrices<Matrix>& matrices,
+ const std::vector<Vector>& u_initial,
+ const std::vector<Vector>& v_initial,
+ const std::vector<Vector>& a_initial);
diff --git a/src/time-stepping/state.cc b/src/time-stepping/state.cc
index 49567d3f..f07a8033 100644
--- a/src/time-stepping/state.cc
+++ b/src/time-stepping/state.cc
@@ -7,9 +7,8 @@
#include "state/sliplawstateupdater.cc"
template <class ScalarVector, class Vector>
-std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(
- Config::stateModel model, ScalarVector const &alpha_initial,
- Dune::BitSetVector<1> const &frictionalNodes, double L, double V0) {
+std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
+ const std::vector<Dune::BitSetVector<_Tp1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0) {
switch (model) {
case Config::AgeingLaw:
return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
diff --git a/src/time-stepping/state.hh b/src/time-stepping/state.hh
index cb372217..18194643 100644
--- a/src/time-stepping/state.hh
+++ b/src/time-stepping/state.hh
@@ -12,7 +12,7 @@
template <class ScalarVector, class Vector>
std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(
- Config::stateModel model, ScalarVector const &alpha_initial,
- Dune::BitSetVector<1> const &frictionalNodes, double L, double V0);
+ Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
+ const std::vector<Dune::BitSetVector<1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);
#endif
diff --git a/src/time-stepping/state/ageinglawstateupdater.cc b/src/time-stepping/state/ageinglawstateupdater.cc
index bbec51f2..24985bc2 100644
--- a/src/time-stepping/state/ageinglawstateupdater.cc
+++ b/src/time-stepping/state/ageinglawstateupdater.cc
@@ -5,8 +5,10 @@
template <class ScalarVector, class Vector>
AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater(
- ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes,
- double _L, double _V0)
+ const std::vector<ScalarVector>& _alpha_initial,
+ const std::vector<Dune::BitSetVector<1>>& _nodes,
+ const std::vector<double>& _L,
+ const std::vector<double>& _V0)
: alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
template <class ScalarVector, class Vector>
@@ -32,21 +34,28 @@ double liftSingularity(double c, double x) {
template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::solve(
- Vector const &velocity_field) {
- for (size_t i = 0; i < nodes.size(); ++i) {
- if (not toBool(nodes[i]))
- continue;
-
- double const V = velocity_field[i].two_norm();
- double const mtoL = -tau / L;
- alpha[i] = std::log(std::exp(alpha_o[i] + V * mtoL) +
- V0 * liftSingularity(mtoL, V));
+ const std::vector<Vector>& velocity_field) {
+ for (size_t i = 0; i < alpha.size(); ++i) {
+ const auto& velocity_fieldi = velocity_field[i];
+ const auto& nodesi = nodes[i];
+ auto& alphai = alpha[i];
+ auto& alpha_oi = alpha_o[i];
+
+ for (size_t j = 0; j < nodesi.size(); ++j) {
+ if (not toBool(nodesi[j]))
+ continue;
+
+ double const V = velocity_fieldi[j].two_norm();
+ double const mtoL = -tau / L[i];
+ alphai[j] = std::log(std::exp(alpha_oi[j] + V * mtoL) +
+ V0[i] * liftSingularity(mtoL, V));
+ }
}
}
template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
- ScalarVector &_alpha) {
+ std::vector<ScalarVector>& _alpha) {
_alpha = alpha;
}
diff --git a/src/time-stepping/state/ageinglawstateupdater.hh b/src/time-stepping/state/ageinglawstateupdater.hh
index 874e26d5..4b8ec558 100644
--- a/src/time-stepping/state/ageinglawstateupdater.hh
+++ b/src/time-stepping/state/ageinglawstateupdater.hh
@@ -6,23 +6,23 @@
template <class ScalarVector, class Vector>
class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
public:
- AgeingLawStateUpdater(ScalarVector const &_alpha_initial,
- Dune::BitSetVector<1> const &_nodes, double _L,
- double _V0);
+ AgeingLawStateUpdater(const std::vector<ScalarVector>& _alpha_initial,
+ const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
+ const std::vector<double>& _V0);
void nextTimeStep() override;
void setup(double _tau) override;
- void solve(Vector const &velocity_field) override;
- void extractAlpha(ScalarVector &) override;
+ void solve(const std::vector<Vector>& velocity_field) override;
+ void extractAlpha(std::vector<ScalarVector> &) override;
std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
private:
- ScalarVector alpha_o;
- ScalarVector alpha;
- Dune::BitSetVector<1> const &nodes;
- double const L;
- double const V0;
+ std::vector<ScalarVector> alpha_o;
+ std::vector<ScalarVector> alpha;
+ const std::vector<Dune::BitSetVector<1>>& nodes;
+ const std::vector<double>& L;
+ const std::vector<double>& V0;
double tau;
};
#endif
diff --git a/src/time-stepping/state/sliplawstateupdater.cc b/src/time-stepping/state/sliplawstateupdater.cc
index a31964d0..6b322d89 100644
--- a/src/time-stepping/state/sliplawstateupdater.cc
+++ b/src/time-stepping/state/sliplawstateupdater.cc
@@ -5,8 +5,9 @@
template <class ScalarVector, class Vector>
SlipLawStateUpdater<ScalarVector, Vector>::SlipLawStateUpdater(
- ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes,
- double _L, double _V0)
+ const std::vector<ScalarVector>& _alpha_initial,
+ const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
+ const std::vector<double>& _V0)
: alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
template <class ScalarVector, class Vector>
@@ -21,21 +22,29 @@ void SlipLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
template <class ScalarVector, class Vector>
void SlipLawStateUpdater<ScalarVector, Vector>::solve(
- Vector const &velocity_field) {
- for (size_t i = 0; i < nodes.size(); ++i) {
- if (not toBool(nodes[i]))
- continue;
-
- double const V = velocity_field[i].two_norm();
- double const mtVoL = -tau * V / L;
- alpha[i] = (V <= 0) ? alpha_o[i] : std::expm1(mtVoL) * std::log(V / V0) +
- alpha_o[i] * std::exp(mtVoL);
- }
+ const std::vector<Vector>& velocity_field) {
+
+ for (size_t i = 0; i < alpha.size(); ++i) {
+ const auto& velocity_fieldi = velocity_field[i];
+ const auto& nodesi = nodes[i];
+ auto& alphai = alpha[i];
+ auto& alpha_oi = alpha_o[i];
+
+ for (size_t j = 0; j < nodesi.size(); ++j) {
+ if (not toBool(nodesi[j]))
+ continue;
+
+ double const V = velocity_fieldi[j].two_norm();
+ double const mtoL = -tau * V / L[i];
+ alphai[j] = (V <= 0) ? alpha_oi[j] : std::expm1(mtVoL) * std::log(V / V0[i]) +
+ alpha_oi[j] * std::exp(mtVoL);
+ }
+ }
}
template <class ScalarVector, class Vector>
void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha(
- ScalarVector &_alpha) {
+ std::vector<ScalarVector>& _alpha) {
_alpha = alpha;
}
diff --git a/src/time-stepping/state/sliplawstateupdater.hh b/src/time-stepping/state/sliplawstateupdater.hh
index ddf2b26e..760044e8 100644
--- a/src/time-stepping/state/sliplawstateupdater.hh
+++ b/src/time-stepping/state/sliplawstateupdater.hh
@@ -6,24 +6,24 @@
template <class ScalarVector, class Vector>
class SlipLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
public:
- SlipLawStateUpdater(ScalarVector const &_alpha_initial,
- Dune::BitSetVector<1> const &_nodes, double _L,
- double _V0);
+ SlipLawStateUpdater(const std::vector<ScalarVector>& _alpha_initial,
+ const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
+ const std::vector<double>& _V0);
- void nextTimeStep() override;
- void setup(double _tau) override;
- void solve(Vector const &velocity_field) override;
- void extractAlpha(ScalarVector &) override;
+ void nextTimeStep() override;
+ void setup(double _tau) override;
+ void solve(const std::vector<Vector>& velocity_field) override;
+ void extractAlpha(std::vector<ScalarVector> &) override;
- std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
+ std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
private:
- ScalarVector alpha_o;
- ScalarVector alpha;
- Dune::BitSetVector<1> const &nodes;
- double const L;
- double const V0;
- double tau;
+ std::vector<ScalarVector> alpha_o;
+ std::vector<ScalarVector> alpha;
+ const std::vector<Dune::BitSetVector<1>>& nodes;
+ const std::vector<double>& L;
+ const std::vector<double>& V0;
+ double tau;
};
#endif
diff --git a/src/time-stepping/state/stateupdater.hh b/src/time-stepping/state/stateupdater.hh
index e1431ac1..f7dd1b1b 100644
--- a/src/time-stepping/state/stateupdater.hh
+++ b/src/time-stepping/state/stateupdater.hh
@@ -7,8 +7,8 @@ template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
void virtual nextTimeStep() = 0;
void virtual setup(double _tau) = 0;
- void virtual solve(Vector const &velocity_field) = 0;
- void virtual extractAlpha(ScalarVector &alpha) = 0;
+ void virtual solve(const std::vector<Vector>& velocity_field) = 0;
+ void virtual extractAlpha(std::vector<ScalarVector>& alpha) = 0;
std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const = 0;
};
diff --git a/src/time-stepping/state_tmpl.cc b/src/time-stepping/state_tmpl.cc
index 20f4b755..5703f81a 100644
--- a/src/time-stepping/state_tmpl.cc
+++ b/src/time-stepping/state_tmpl.cc
@@ -2,5 +2,5 @@
template std::shared_ptr<StateUpdater<ScalarVector, Vector>>
initStateUpdater<ScalarVector, Vector>(
- Config::stateModel model, ScalarVector const &alpha_initial,
- Dune::BitSetVector<1> const &frictionalNodes, double L, double V0);
+ Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
+ const std::vector<Dune::BitSetVector<_Tp1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);
--
GitLab