diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh
index 9378c10677101b25e94f213c216b9fcdf41f1fed..5ee5bc9ed85b42eccdda19a2569d11af8d03c37f 100644
--- a/dune/tectonic/globalruinanonlinearity.hh
+++ b/dune/tectonic/globalruinanonlinearity.hh
@@ -26,13 +26,14 @@ class GlobalRuinaNonlinearity
shared_ptr<BlockVector<FieldVector<double, 1>> const> mu,
shared_ptr<BlockVector<FieldVector<double, 1>> const> eta,
shared_ptr<BlockVector<FieldVector<double, 1>> const> normalStress,
- shared_ptr<BlockVector<FieldVector<double, 1>> const> state)
+ shared_ptr<BlockVector<FieldVector<double, 1>> const> state, double h)
: nodalIntegrals(nodalIntegrals),
a(a),
mu(mu),
eta(eta),
normalStress(normalStress),
state(state),
+ h(h),
trivialNonlinearity(
new LocalNonlinearity<dim>(make_shared<TrivialFunction const>())),
restrictions(nodalIntegrals->size()) // TODO: can we get the size from
@@ -54,7 +55,7 @@ class GlobalRuinaNonlinearity
auto const func = make_shared<RuinaFunction const>(
(*nodalIntegrals)[i][0], (*a)[i][0], (*mu)[i][0], (*eta)[i][0],
- (*normalStress)[i][0], (*state)[i][0]);
+ (*normalStress)[i][0], (*state)[i][0], h);
restrictions[i] = make_shared<LocalNonlinearity<dim> const>(func);
return restrictions[i];
}
@@ -68,6 +69,7 @@ class GlobalRuinaNonlinearity
shared_ptr<BlockVector<FieldVector<double, 1>> const> eta;
shared_ptr<BlockVector<FieldVector<double, 1>> const> normalStress;
shared_ptr<BlockVector<FieldVector<double, 1>> const> state;
+ double const h;
std::vector<shared_ptr<LocalNonlinearity<dim> const>> mutable restrictions;
};
diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/nicefunction.hh
index d144d3ec7725a6d59a8768ac151ff01c6f7cabc7..015bb5025e61fadb54e02c761a442ba4a3c87dec 100644
--- a/dune/tectonic/nicefunction.hh
+++ b/dune/tectonic/nicefunction.hh
@@ -31,14 +31,15 @@ class NiceFunction : public VirtualFunction<double, double> {
class RuinaFunction : public NiceFunction {
public:
RuinaFunction(double coefficient, double a, double mu, double eta,
- double normalStress, double state)
+ double normalStress, double state, double h)
: coefficient(coefficient),
a(a),
mu(mu),
rho(exp(-mu / a)),
eta(eta),
normalStress(normalStress),
- state(state) {}
+ state(state),
+ h(h) {}
/*
If mu and sigma_n denote the coefficient of friction and the
@@ -48,17 +49,18 @@ class RuinaFunction : public NiceFunction {
with the constants a and b from Ruina's model and
- h(beta) = beta (a (log beta - 1) + mu)
+ r(beta) = beta (a (log beta - 1) + mu)
as well as
rho = exp(-mu/a)
*/
void virtual evaluate(double const &x, double &y) const {
- double const arg = std::max(eta * x, rho);
- double const h = arg * (a * (std::log(arg) - 1) + mu + state);
- y = 1 / eta * normalStress * h + a * rho + mu + state;
+ double const arg = std::max(eta * x / h, rho);
+ double const r = arg * (a * (std::log(arg) - 1) + mu + state);
+ y = 1 / eta * normalStress * r + a * rho + mu + state;
y *= coefficient;
+ y *= h;
}
/*
@@ -69,18 +71,20 @@ class RuinaFunction : public NiceFunction {
= a * log(eta x) + mu
*/
double virtual leftDifferential(double s) const {
- if (eta * s <= rho)
+ double const arg = eta * s / h;
+ if (arg <= rho)
return 0;
- return coefficient * normalStress * (a * std::log(eta * s) + mu + state);
+ return coefficient * normalStress * (a * std::log(arg) + mu + state);
}
/* see above */
double virtual rightDifferential(double s) const {
- if (eta * s <= rho)
+ double const arg = eta * s / h;
+ if (arg <= rho)
return 0;
- return coefficient * normalStress * (a * std::log(eta * s) + mu + state);
+ return coefficient * normalStress * (a * std::log(arg) + mu + state);
}
/*
@@ -90,10 +94,10 @@ class RuinaFunction : public NiceFunction {
*/
double virtual second_deriv(double s) const {
// includes the case eta * s = rho for which there is no second derivative
- if (eta * s <= rho)
+ if (eta * s / h <= rho)
return 0;
else
- return coefficient * normalStress * (a / s);
+ return coefficient * normalStress * a / s;
}
double virtual regularity(double s) const {
@@ -110,6 +114,7 @@ class RuinaFunction : public NiceFunction {
double const eta;
double const normalStress;
double const state;
+ double const h;
double const rho;
};
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 95c0f05e26b4745bd443ee6215f5e5f87876d564..c403b30bfa7d612cfda82bf95267ac19f2e36b90 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -103,11 +103,12 @@ template <class GridType, class GridView, class LocalVectorType, class FEBasis>
void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
Dune::BitSetVector<1> const &neumannNodes,
Dune::BlockVector<LocalVectorType> &f,
- int run) { // constant sample function on neumann boundary
+ double time) { // constant sample function on neumann
+ // boundary
BoundaryPatch<GridView> neumannBoundary(gridView, neumannNodes);
LocalVectorType SampleVector(0);
// FIXME: random values (time-dependent)
- SampleVector[0] = run / 10.0;
+ SampleVector[0] = time;
SampleVector[1] = 0;
ConstantFunction<LocalVectorType, LocalVectorType> fNeumann(SampleVector);
NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
@@ -145,7 +146,8 @@ assemble_nonlinearity(
int size, Dune::ParameterTree const &parset,
Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
nodalIntegrals,
- Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state) {
+ Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
+ double h) {
typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
// {{{ Assemble terms for the nonlinearity
@@ -166,7 +168,7 @@ assemble_nonlinearity(
return Dune::make_shared<
Dune::GlobalRuinaNonlinearity<VectorType, MatrixType> const>(
- nodalIntegrals, a, mu, eta, normalStress, state);
+ nodalIntegrals, a, mu, eta, normalStress, state, h);
} else if (friction_model == std::string("Laursen")) {
return Dune::make_shared<Dune::GlobalLaursenNonlinearity<
Dune::LinearFunction, VectorType, MatrixType> const>(mu, normalStress,
@@ -312,6 +314,7 @@ int main(int argc, char *argv[]) {
timer.reset();
auto const timesteps = parset.get<size_t>("timesteps");
+ double h = 1.0 / timesteps;
for (size_t run = 1; run <= timesteps; ++run) {
if (parset.get<bool>("printProgress")) {
std::cout << ".";
@@ -320,7 +323,7 @@ int main(int argc, char *argv[]) {
// b = neumann
assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
- leafView, p1Basis, neumannNodes, b1, run);
+ leafView, p1Basis, neumannNodes, b1, h * run);
b2 = b1;
b3 = b1;
b4 = b1;
@@ -339,7 +342,7 @@ int main(int argc, char *argv[]) {
auto myGlobalNonlinearity =
assemble_nonlinearity<VectorType, OperatorType>(
grid->size(grid->maxLevel(), dim), parset, nodalIntegrals,
- state);
+ state, h);
MyConvexProblemType myConvexProblem(stiffnessMatrix,
*myGlobalNonlinearity, b1);
MyBlockProblemType myBlockProblem(parset, myConvexProblem);
@@ -362,7 +365,7 @@ int main(int argc, char *argv[]) {
auto myGlobalNonlinearity =
assemble_nonlinearity<VectorType, OperatorType>(
grid->size(grid->maxLevel(), dim), parset, nodalIntegrals,
- state);
+ state, h);
MyConvexProblemType myConvexProblem(stiffnessMatrix,
*myGlobalNonlinearity, b4);