diff --git a/dune/tectonic/frictiondata.hh b/dune/tectonic/frictiondata.hh
index 0262f2df8529be50b73e3dd8d1a3caf89596062b..5f2c43c66a2554ee619f7e32c290b20967c83a04 100644
--- a/dune/tectonic/frictiondata.hh
+++ b/dune/tectonic/frictiondata.hh
@@ -3,12 +3,8 @@
#include <dune/common/parametertree.hh>
struct FrictionData {
- FrictionData(Dune::ParameterTree const &parset)
- : L(parset.get<double>("L")),
- V0(parset.get<double>("V0")),
- a(parset.get<double>("a")),
- b(parset.get<double>("b")),
- mu0(parset.get<double>("mu0")) {}
+ explicit FrictionData(double _L, double _V0, double _a, double _b, double _mu0)
+ : L(_L), V0(_V0), a(_a), b(_b), mu0(_mu0) {}
double const L;
double const V0;
diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index eb447bfb4880560633e469fbda7c2b430aae83b1..03d74518ddb75ca0afd1994f172f8a806f7b8983 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -18,6 +18,7 @@ class FrictionPotentialWrapper {
double virtual differential(double s) const = 0;
double virtual second_deriv(double x) const = 0;
double virtual regularity(double s) const = 0;
+ double virtual coefficientOfFriction(double s) const = 0;
double virtual evaluate(double x) const {
DUNE_THROW(Dune::NotImplemented, "evaluation not implemented");
@@ -28,16 +29,19 @@ class FrictionPotentialWrapper {
class FrictionPotential : public FrictionPotentialWrapper {
public:
- FrictionPotential(double coefficient, double _normalStress,
- FrictionData const &_fd)
+ FrictionPotential(double coefficient, double _normalStress, FrictionData _fd)
: fd(_fd), weight(coefficient), normalStress(_normalStress) {}
- double differential(double V) const {
+ double coefficientOfFriction(double V) const {
assert(V >= 0.0);
if (V <= V_cutoff)
return 0.0;
- return weight * (-normalStress) * fd.a * (std::log(V) - logV_m);
+ return fd.a * (std::log(V) - logV_m);
+ }
+
+ double differential(double V) const {
+ return weight * (-normalStress) * coefficientOfFriction(V);
}
double second_deriv(double V) const {
@@ -64,7 +68,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
}
private:
- FrictionData const &fd;
+ FrictionData const fd;
double const weight;
double const normalStress;
double logV_m;
@@ -75,6 +79,8 @@ class TrivialFunction : public FrictionPotentialWrapper {
public:
double evaluate(double) const { return 0; }
+ double coefficientOfFriction(double s) const { return 0; }
+
double differential(double) const { return 0; }
double second_deriv(double) const { return 0; }
diff --git a/dune/tectonic/globalfrictiondata.hh b/dune/tectonic/globalfrictiondata.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6500352399d898aa8139c76d4b88bea7548ab665
--- /dev/null
+++ b/dune/tectonic/globalfrictiondata.hh
@@ -0,0 +1,36 @@
+#ifndef GLOBAL_FRICTION_DATA_HH
+#define GLOBAL_FRICTION_DATA_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/tectonic/frictiondata.hh>
+
+template <class DomainType>
+double evaluateScalarFunction(
+ Dune::VirtualFunction<DomainType, Dune::FieldVector<double, 1>> const &f,
+ DomainType const &x) {
+ Dune::FieldVector<double, 1> ret;
+ f.evaluate(x, ret);
+ return ret;
+};
+
+template <int dimension> class GlobalFrictionData {
+public:
+ FrictionData operator()(Dune::FieldVector<double, dimension> const &x) const {
+ return FrictionData(L(), V0(), evaluateScalarFunction(a(), x),
+ evaluateScalarFunction(b(), x), mu0());
+ }
+
+protected:
+ using VirtualFunction = Dune::VirtualFunction<
+ Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
+
+ double virtual const &L() const = 0;
+ double virtual const &V0() const = 0;
+ VirtualFunction virtual const &a() const = 0;
+ VirtualFunction virtual const &b() const = 0;
+ double virtual const &mu0() const = 0;
+};
+
+#endif
diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh
index c95ed1f7a2bd407cc725efb4771703055cecb3be..4dc996d805c8aa7d9c78dbe7c5faf81728bc9ed3 100644
--- a/dune/tectonic/globalruinanonlinearity.hh
+++ b/dune/tectonic/globalruinanonlinearity.hh
@@ -6,13 +6,16 @@
#include <dune/common/bitsetvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
+#include <dune/grid/common/mcmgmapper.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
#include "globalnonlinearity.hh"
#include "frictionpotential.hh"
-template <class Matrix, class Vector>
+template <class Matrix, class Vector, class GridView>
class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
public:
using GlobalNonlinearity<Matrix, Vector>::block_size;
@@ -22,19 +25,28 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
using typename GlobalNonlinearity<Matrix, Vector>::ScalarVector;
public:
- GlobalRuinaNonlinearity(Dune::BitSetVector<1> const &frictionalNodes,
- ScalarVector const &weights, FrictionData const &fd,
+ GlobalRuinaNonlinearity(BoundaryPatch<GridView> const &frictionalBoundary,
+ GridView const &gridView,
+ GlobalFrictionData<block_size> const &frictionInfo,
+ // Note: passing the following two makes no sense
+ ScalarVector const &weights,
ScalarVector const &normalStress)
- : restrictions(frictionalNodes.size()) {
+ : restrictions(normalStress.size()) {
auto trivialNonlinearity =
std::make_shared<Friction>(std::make_shared<TrivialFunction>());
- for (size_t i = 0; i < restrictions.size(); ++i) {
- if (not frictionalNodes[i][0]) {
+
+ Dune::MultipleCodimMultipleGeomTypeMapper<
+ GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
+ for (auto it = gridView.template begin<block_size>();
+ it != gridView.template end<block_size>(); ++it) {
+ auto const i = vertexMapper.map(*it);
+ auto const coordinate = it->geometry().corner(0);
+ if (not frictionalBoundary.containsVertex(i)) {
restrictions[i] = trivialNonlinearity;
continue;
}
- auto const fp =
- std::make_shared<FrictionPotential>(weights[i], normalStress[i], fd);
+ auto const fp = std::make_shared<FrictionPotential>(
+ weights[i], normalStress[i], frictionInfo(coordinate));
restrictions[i] = std::make_shared<Friction>(fp);
}
}
@@ -44,6 +56,13 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
restrictions[i]->updateLogState(logState[i]);
}
+ void coefficientOfFriction(Vector const &x, ScalarVector &coefficient) {
+ assert(x.size() == restrictions.size());
+ coefficient.resize(restrictions.size());
+ for (size_t i = 0; i < restrictions.size(); ++i)
+ coefficient[i] = restrictions[i]->coefficientOfFriction(x[i]);
+ }
+
/*
Return a restriction of the outer function to the i'th node.
*/
diff --git a/dune/tectonic/localfriction.hh b/dune/tectonic/localfriction.hh
index 04b8121a9b3d06c82a15c44bed2fa4f2471ef67c..5c67dc38e1b496ce2f0306c3f8fb1e73521b0cbd 100644
--- a/dune/tectonic/localfriction.hh
+++ b/dune/tectonic/localfriction.hh
@@ -51,6 +51,10 @@ template <size_t dimension> class LocalFriction {
return func->regularity(xnorm);
}
+ double coefficientOfFriction(VectorType const &x) const {
+ return func->coefficientOfFriction(x.two_norm());
+ }
+
// directional subdifferential: at u on the line u + t*v
// u and v are assumed to be non-zero
void directionalSubDiff(VectorType const &x, VectorType const &v,
diff --git a/src/assemblers.cc b/src/assemblers.cc
index 205ef090a512c19cd060b7eead81ff729970a0f7..b9a3e76b93f0ebf96e011b61cbbda9b36d6b847a 100644
--- a/src/assemblers.cc
+++ b/src/assemblers.cc
@@ -92,9 +92,10 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
template <class GridView, int dimension>
auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
- BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd,
+ BoundaryPatch<GridView> const &frictionalBoundary,
+ GlobalFrictionData<dimension> const &frictionInfo,
ScalarVector const &normalStress)
- -> std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>> {
+ -> std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector, GridView>> {
// Lump negative normal stress (kludge)
ScalarVector weights;
{
@@ -105,13 +106,8 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
weights, frictionalBoundary);
}
- for (size_t i = 0; i < weights.size(); ++i)
- assert(weights[i] >= 0.0);
-
- Dune::BitSetVector<1> frictionalNodes;
- frictionalBoundary.getVertices(frictionalNodes);
- return std::make_shared<GlobalRuinaNonlinearity<Matrix, Vector>>(
- frictionalNodes, weights, fd, normalStress);
+ return std::make_shared<GlobalRuinaNonlinearity<Matrix, Vector, GridView>>(
+ frictionalBoundary, gridView, frictionInfo, weights, normalStress);
}
template <class GridView, int dimension>
diff --git a/src/assemblers.hh b/src/assemblers.hh
index 0672350058e848a6a766518a83c9356c8bf6973d..a281251e7f91d76c08c12734253bd2b5a876e13b 100644
--- a/src/assemblers.hh
+++ b/src/assemblers.hh
@@ -13,6 +13,7 @@
#pragma clang diagnostic pop
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
#include <dune/tectonic/globalruinanonlinearity.hh>
template <class GridView, int dimension> class MyAssembler {
@@ -66,9 +67,10 @@ template <class GridView, int dimension> class MyAssembler {
Dune::VirtualFunction<double, double> const &neumann,
double relativeTime);
- std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>>
+ std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector, GridView>>
assembleFrictionNonlinearity(
- BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd,
+ BoundaryPatch<GridView> const &frictionalBoundary,
+ GlobalFrictionData<dimension> const &frictionInfo,
ScalarVector const &normalStress);
void assembleVonMisesStress(double youngModulus, double poissonRatio,
diff --git a/src/friction_writer.cc b/src/friction_writer.cc
index 46ee0e5bee3c8411c470cc4721484c983282ff6d..22cfe296f887343a1fb0121cd111e1d02ad51455 100644
--- a/src/friction_writer.cc
+++ b/src/friction_writer.cc
@@ -4,20 +4,12 @@
#include "friction_writer.hh"
-double computeCOF(FrictionData const &fd, double V, double logState) {
- double const mu = fd.mu0 + fd.a * std::log(V / fd.V0) +
- fd.b * (logState + std::log(fd.V0 / fd.L));
- return (mu > 0) ? mu : 0;
-}
-
template <class BitVector>
-FrictionWriter<BitVector>::FrictionWriter(FrictionData const &_fd,
- BitVector const &_boundaryNodes)
+FrictionWriter<BitVector>::FrictionWriter(BitVector const &_boundaryNodes)
: coefficientWriter("coefficients", std::fstream::out),
displacementWriter("displacements", std::fstream::out),
stateWriter("states", std::fstream::out),
velocityWriter("velocities", std::fstream::out),
- fd(_fd),
boundaryNodes(_boundaryNodes) {}
template <class BitVector> FrictionWriter<BitVector>::~FrictionWriter() {
@@ -29,13 +21,14 @@ template <class BitVector> FrictionWriter<BitVector>::~FrictionWriter() {
template <class BitVector>
template <class ScalarVector, class Vector>
-void FrictionWriter<BitVector>::writeInfo(ScalarVector const &alpha,
+void FrictionWriter<BitVector>::writeInfo(ScalarVector const &coefficient,
+ ScalarVector const &alpha,
Vector const &u, Vector const &v) {
for (size_t i = 0; i < boundaryNodes.size(); ++i) {
if (!boundaryNodes[i][0])
continue;
- coefficientWriter << computeCOF(fd, v[i].two_norm(), alpha[i]) << " ";
+ coefficientWriter << coefficient[i] << " ";
displacementWriter << u[i][0] << " ";
stateWriter << alpha[i][0] << " ";
velocityWriter << v[i][0] << " ";
diff --git a/src/friction_writer.hh b/src/friction_writer.hh
index c2af5fffa9a7f893ff113c7b902ed776781c2b48..7494b890537a9eb8f83af0828fdec2ee4ea386d8 100644
--- a/src/friction_writer.hh
+++ b/src/friction_writer.hh
@@ -3,23 +3,21 @@
#include <fstream>
-#include <dune/tectonic/frictiondata.hh>
-
template <class BitVector> class FrictionWriter {
public:
- FrictionWriter(FrictionData const &_fd, BitVector const &_boundaryNodes);
+ FrictionWriter(BitVector const &_boundaryNodes);
~FrictionWriter();
template <class ScalarVector, class Vector>
- void writeInfo(ScalarVector const &alpha, Vector const &u, Vector const &v);
+ void writeInfo(ScalarVector const &coefficient, ScalarVector const &alpha,
+ Vector const &u, Vector const &v);
private:
std::fstream coefficientWriter;
std::fstream displacementWriter;
std::fstream stateWriter;
std::fstream velocityWriter;
- FrictionData const &fd;
BitVector const &boundaryNodes;
};
#endif
diff --git a/src/friction_writer_tmpl.cc b/src/friction_writer_tmpl.cc
index dbeb7b8b8c5cabb9eda1639f0cd96ab4468aa34c..9a95b2611f6977c56650d6018cfa92514caa1c00 100644
--- a/src/friction_writer_tmpl.cc
+++ b/src/friction_writer_tmpl.cc
@@ -10,6 +10,6 @@ using BitVector = Dune::BitSetVector<1>;
template class FrictionWriter<BitVector>;
-template void FrictionWriter<BitVector>::writeInfo(ScalarVector const &alpha,
- Vector const &u,
- Vector const &v);
+template void FrictionWriter<BitVector>::writeInfo(
+ ScalarVector const &coefficient, ScalarVector const &alpha, Vector const &u,
+ Vector const &v);
diff --git a/src/mybody.hh b/src/mybody.hh
index ec735642edbca78cf2e0d63f02421c2f4131da77..3b33c93fa023bf208a4bde66327b75d3fa2d49d8 100644
--- a/src/mybody.hh
+++ b/src/mybody.hh
@@ -13,7 +13,6 @@
template <int dimension> class MyBody : public Body<dimension> {
using typename Body<dimension>::ScalarFunction;
using typename Body<dimension>::VectorField;
-
using MyConstantFunction = ConstantFunction<
Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
diff --git a/src/myglobalfrictiondata.hh b/src/myglobalfrictiondata.hh
new file mode 100644
index 0000000000000000000000000000000000000000..43d6b9cfa25d5a697d3f3fc57b0ed8b0a3242515
--- /dev/null
+++ b/src/myglobalfrictiondata.hh
@@ -0,0 +1,43 @@
+#ifndef MY_GLOBAL_FRICTION_DATA_HH
+#define MY_GLOBAL_FRICTION_DATA_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/functions/constantfunction.hh>
+
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "mygeometry.hh"
+
+template <int dimension>
+class MyGlobalFrictionData : public GlobalFrictionData<dimension> {
+private:
+ using typename GlobalFrictionData<dimension>::VirtualFunction;
+
+public:
+ MyGlobalFrictionData(Dune::ParameterTree const &parset, MyGeometry const &tri)
+ : L_(parset.get<double>("L")),
+ V0_(parset.get<double>("V0")),
+ a_(parset.get<double>("a")),
+ b_(parset.get<double>("b")),
+ mu0_(parset.get<double>("mu0")) {}
+
+ double virtual const &L() const override { return L_; }
+ double virtual const &V0() const override { return V0_; }
+ VirtualFunction virtual const &a() const override { return a_; }
+ VirtualFunction virtual const &b() const override { return b_; }
+ double virtual const &mu0() const override { return mu0_; }
+
+private:
+ using MyConstantFunction = ConstantFunction<
+ Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
+
+ double const L_;
+ double const V0_;
+ MyConstantFunction const a_;
+ MyConstantFunction const b_;
+ double const mu0_;
+};
+
+#endif
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 651553c44ede700d2f321dc1ac2a4cd7e3c08e3e..0a10471c91cdb280817eebd6ad6c3a99c337ad48 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -76,6 +76,7 @@
#include "friction_writer.hh"
#include "mybody.hh"
#include "mygeometry.hh"
+#include "myglobalfrictiondata.hh"
#include "mygrid.hh"
#include "solverfactory.hh"
#include "state.hh"
@@ -193,14 +194,15 @@ int main(int argc, char *argv[]) {
normalStress = (myGeometry.A[1] - myGeometry.C[1]) *
parset.get<double>("gravity") *
parset.get<double>("body.density");
- FrictionData const frictionData(parset.sub("boundary.friction"));
+ MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"),
+ myGeometry);
// {{{ Initial conditions
ScalarVector alpha_initial(fineVertexCount);
alpha_initial =
std::log(parset.get<double>("boundary.friction.initialState"));
auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
- frictionalBoundary, frictionData, normalStress);
+ frictionalBoundary, frictionInfo, normalStress);
myGlobalNonlinearity->updateLogState(alpha_initial);
using LinearFactory = SolverFactory<
@@ -267,8 +269,12 @@ int main(int argc, char *argv[]) {
parset.sub("a0.solver"));
}
// }}}
- FrictionWriter<Dune::BitSetVector<1>> writer(frictionData, frictionalNodes);
- writer.writeInfo(alpha_initial, u_initial, v_initial);
+ FrictionWriter<Dune::BitSetVector<1>> writer(frictionalNodes);
+ {
+ ScalarVector c;
+ myGlobalNonlinearity->coefficientOfFriction(v_initial, c);
+ writer.writeInfo(c, alpha_initial, u_initial, v_initial);
+ }
MyVTKWriter<typename MyAssembler::VertexBasis,
typename MyAssembler::CellBasis> const
@@ -310,7 +316,8 @@ int main(int argc, char *argv[]) {
u_initial, v_initial, a_initial);
auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
- alpha_initial, frictionalNodes, frictionData);
+ alpha_initial, frictionalNodes,
+ parset.get<double>("boundary.friction.L"));
Vector v = v_initial;
ScalarVector alpha(fineVertexCount);
@@ -418,7 +425,11 @@ int main(int argc, char *argv[]) {
if (printProgress)
std::cout << std::endl;
- writer.writeInfo(alpha, u, v);
+ {
+ ScalarVector c;
+ myGlobalNonlinearity->coefficientOfFriction(v, c);
+ writer.writeInfo(c, alpha, u, v);
+ }
iterationWriter << std::endl;
relaxationWriter << std::endl;
diff --git a/src/state.hh b/src/state.hh
index e1d840201fd226fbf8c01caaf3c88572c2f2c139..779640d0a0ff102441e4c653c4dd28a38fd57339 100644
--- a/src/state.hh
+++ b/src/state.hh
@@ -11,14 +11,14 @@
template <class ScalarVector, class Vector>
std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(
Config::stateModel model, ScalarVector const &alpha_initial,
- Dune::BitSetVector<1> const &frictionalNodes, FrictionData const &fd) {
+ Dune::BitSetVector<1> const &frictionalNodes, double L) {
switch (model) {
case Config::Dieterich:
return std::make_shared<DieterichStateUpdater<ScalarVector, Vector>>(
- alpha_initial, frictionalNodes, fd.L);
+ alpha_initial, frictionalNodes, L);
case Config::Ruina:
return std::make_shared<RuinaStateUpdater<ScalarVector, Vector>>(
- alpha_initial, frictionalNodes, fd.L);
+ alpha_initial, frictionalNodes, L);
default:
assert(false);
}