diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh
index 2cac0eb69c19147be3ee63367799dd7f9734671f..257b472f8785cc97e92f1fc76eaf998d4c6dfa87 100644
--- a/dune/tectonic/globalruinanonlinearity.hh
+++ b/dune/tectonic/globalruinanonlinearity.hh
@@ -25,18 +25,17 @@ class GlobalRuinaNonlinearity
GlobalRuinaNonlinearity(dataref nodalIntegrals, dataref a, dataref mu,
dataref eta, dataref normalStress, dataref b,
- dataref state, dataref L, double h)
+ dataref state, dataref L)
: restrictions(nodalIntegrals.size()) {
auto trivialNonlinearity = make_shared<LocalNonlinearity<dim> const>(
make_shared<TrivialFunction const>());
for (size_t i = 0; i < restrictions.size(); ++i) {
- restrictions[i] =
- nodalIntegrals[i] == 0
- ? trivialNonlinearity
- : make_shared<LocalNonlinearity<dim> const>(
- make_shared<RuinaFunction const>(
- nodalIntegrals[i], a[i], mu[i], eta[i], normalStress[i],
- b[i], state[i], L[i], h));
+ restrictions[i] = nodalIntegrals[i] == 0
+ ? trivialNonlinearity
+ : make_shared<LocalNonlinearity<dim> const>(
+ make_shared<RuinaFunction const>(
+ nodalIntegrals[i], a[i], mu[i], eta[i],
+ normalStress[i], b[i], state[i], L[i]));
}
}
diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/nicefunction.hh
index 10b7417d64510756337e924ac099eb6a9a99b8b6..f1647d5a7bc24a25693601253c360223e6041322 100644
--- a/dune/tectonic/nicefunction.hh
+++ b/dune/tectonic/nicefunction.hh
@@ -41,18 +41,17 @@ class NiceFunction {
class RuinaFunction : public NiceFunction {
public:
RuinaFunction(double coefficient, double a, double mu, double eta,
- double normalStress, double b, double state, double L, double h)
+ double normalStress, double b, double state, double L)
: NiceFunction(),
a(a),
eta(eta),
- h(h),
coefficientProduct(coefficient * normalStress),
K(mu + b * (state - std::log(eta * L))), // state is assumed to be
// logarithmic
rho(std::exp(-K / a)) {}
double virtual evaluate(double x) const {
- double const arg = x / h * eta;
+ double const arg = x * eta;
if (arg <= rho)
return 0;
@@ -61,7 +60,7 @@ class RuinaFunction : public NiceFunction {
}
double virtual leftDifferential(double x) const {
- double const arg = x / h * eta;
+ double const arg = x * eta;
if (arg <= rho)
return 0;
@@ -74,7 +73,7 @@ class RuinaFunction : public NiceFunction {
}
double virtual second_deriv(double x) const {
- double const arg = x / h * eta;
+ double const arg = x * eta;
if (arg <= rho)
return 0;
@@ -82,7 +81,7 @@ class RuinaFunction : public NiceFunction {
}
double virtual regularity(double x) const {
- double const arg = x / h * eta;
+ double const arg = x * eta;
// TODO: Make this controllable
if (std::abs(arg - rho) < 1e-14)
return std::numeric_limits<double>::infinity();
@@ -95,7 +94,6 @@ class RuinaFunction : public NiceFunction {
private:
double const a;
double const eta;
- double const h;
double const coefficientProduct;
double const K;
double const rho;
diff --git a/src/assemblers.cc b/src/assemblers.cc
index cf215ea3c81bbb170001c81a6908ab01065cbe7c..b027e72fa91258ad4c06b05d7a1e63a0bf659a95 100644
--- a/src/assemblers.cc
+++ b/src/assemblers.cc
@@ -84,7 +84,7 @@ assemble_nonlinearity(
return Dune::make_shared<
Dune::GlobalRuinaNonlinearity<MatrixType, VectorType> const>(
- nodalIntegrals, a, mu, eta, normalStress, b, state, L, tau);
+ nodalIntegrals, a, mu, eta, normalStress, b, state, L);
}
case Config::Laursen:
assert(false);
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index e638d0e80b6c9c78d34c70157f417b49b06a1c07..1b6ce4d5dd939d461e7701c59c04b5e553467f64 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -270,11 +270,11 @@ int main(int argc, char *argv[]) {
if (run == 1 || !parset.get<bool>("implicitTwoStep"))
timeSteppingScheme =
Dune::make_shared<IE>(ell, stiffnessMatrix, u, u_old_old,
- ignoreNodes, dirichletFunction, time);
+ ignoreNodes, dirichletFunction, time, tau);
else
timeSteppingScheme =
Dune::make_shared<ITS>(ell, stiffnessMatrix, u, u_old_old,
- ignoreNodes, dirichletFunction, time);
+ ignoreNodes, dirichletFunction, time, tau);
timeSteppingScheme->setup(problem_rhs, problem_iterate, problem_A);
diff --git a/src/timestepping.cc b/src/timestepping.cc
index 71d8aae094ac8d25e61c418a0f7821b74e736368..a00e622f879914d42fcb4fe0ddcc7b5cb509241c 100644
--- a/src/timestepping.cc
+++ b/src/timestepping.cc
@@ -4,14 +4,16 @@ class TimeSteppingScheme {
TimeSteppingScheme(VectorType const &_ell, MatrixType const &_A,
VectorType const &_u_old, VectorType const *_u_old_old,
Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction, double _time)
+ FunctionType const &_dirichletFunction, double _time,
+ double _tau)
: ell(_ell),
A(_A),
u_old(_u_old),
u_old_old(_u_old_old),
dirichletNodes(_dirichletNodes),
dirichletFunction(_dirichletFunction),
- time(_time) {}
+ time(_time),
+ tau(_tau) {}
virtual ~TimeSteppingScheme() {}
@@ -31,7 +33,8 @@ class TimeSteppingScheme {
VectorType const *u_old_old;
Dune::BitSetVector<dim> const &dirichletNodes;
FunctionType const &dirichletFunction;
- double time;
+ double const time;
+ double const tau;
};
// Implicit Euler: Solve the problem
@@ -50,10 +53,13 @@ class ImplicitEuler
void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
MatrixType &problem_A) const {
- problem_A = this->A;
problem_rhs = this->ell;
- problem_A.mmv(this->u_old, problem_rhs);
+ this->A.mmv(this->u_old, problem_rhs);
+
+ problem_A = this->A;
+ problem_A *= this->tau;
+ // Treat as Delta u_n for now
problem_iterate = this->u_old;
if (this->u_old_old)
problem_iterate -= *this->u_old_old;
@@ -67,17 +73,21 @@ class ImplicitEuler
val - this->u_old[i][0]; // Time-dependent X direction
} else if (this->dirichletNodes[i][1])
problem_iterate[i][1] = 0; // Y direction described
+
+ // Make it a velocity
+ problem_iterate /= this->tau;
}
void virtual extractSolution(VectorType const &problem_iterate,
VectorType &solution) const {
solution = this->u_old;
- solution += problem_iterate;
+ solution.axpy(this->tau, problem_iterate);
}
void virtual extractDifference(VectorType const &problem_iterate,
VectorType &difference) const {
difference = problem_iterate;
+ difference *= this->tau;
}
};
@@ -94,13 +104,14 @@ class ImplicitTwoStep
void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
MatrixType &problem_A) const {
- problem_A = this->A;
- problem_A /= 1.5;
problem_rhs = this->ell;
- problem_A.usmv(-2, this->u_old, problem_rhs);
- problem_A.usmv(.5, *this->u_old_old, problem_rhs);
+ this->A.usmv(-4.0 / 3.0, this->u_old, problem_rhs);
+ this->A.usmv(+1.0 / 3.0, *this->u_old_old, problem_rhs);
+
+ problem_A = this->A;
+ problem_A *= 2.0 / 3.0 * this->tau;
- // The finite difference makes a good start
+ // The finite difference makes a good start (treat it as such for now)
problem_iterate = this->u_old;
problem_iterate -= *this->u_old_old;
@@ -117,11 +128,14 @@ class ImplicitTwoStep
// Y direction described
problem_iterate[i][1] =
-2 * this->u_old[i][1] + .5 * (*this->u_old_old)[i][1];
+ // Now treat it as a velocity
+ problem_iterate /= this->tau;
}
void virtual extractSolution(VectorType const &problem_iterate,
VectorType &solution) const {
solution = problem_iterate;
+ solution *= this->tau;
solution.axpy(2, this->u_old);
solution.axpy(-.5, *this->u_old_old);
solution *= 2.0 / 3.0;
@@ -129,6 +143,7 @@ class ImplicitTwoStep
// Check if we split correctly
{
VectorType test = problem_iterate;
+ test *= this->tau;
test.axpy(-1.5, solution);
test.axpy(+2, this->u_old);
test.axpy(-.5, *this->u_old_old);
@@ -139,5 +154,6 @@ class ImplicitTwoStep
void virtual extractDifference(VectorType const &problem_iterate,
VectorType &difference) const {
difference = problem_iterate;
+ difference *= this->tau;
}
};