From d49dfbc17af9e81b28489eb79972e60c7a8b636c Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Wed, 14 Dec 2011 14:54:33 +0100
Subject: [PATCH] Handle irregularity
---
dune/tectonic/globalnonlinearity.hh | 5 +++++
dune/tectonic/localnonlinearity.hh | 12 ++++++++++--
dune/tectonic/myblockproblem.hh | 19 +++----------------
dune/tectonic/nicefunction.hh | 20 +++++++++++++++++---
src/one-body-sample.parset | 2 +-
5 files changed, 36 insertions(+), 22 deletions(-)
diff --git a/dune/tectonic/globalnonlinearity.hh b/dune/tectonic/globalnonlinearity.hh
index dc773dfc..9a284eaa 100644
--- a/dune/tectonic/globalnonlinearity.hh
+++ b/dune/tectonic/globalnonlinearity.hh
@@ -33,6 +33,11 @@ class GlobalNonlinearity {
res->addGradient(v[i], gradient[i]);
}
}
+
+ double regularity(int i, typename VectorType::block_type const &x) const {
+ auto res = restriction(i);
+ return res->regularity(x);
+ }
};
}
#endif
diff --git a/dune/tectonic/localnonlinearity.hh b/dune/tectonic/localnonlinearity.hh
index a6ce8e3b..81c81552 100644
--- a/dune/tectonic/localnonlinearity.hh
+++ b/dune/tectonic/localnonlinearity.hh
@@ -26,6 +26,13 @@ template <int dimension> class LocalNonlinearity {
return ret;
}
+ double regularity(VectorType const &x) const {
+ if (x.two_norm() < 1e-8) // TODO
+ return std::numeric_limits<double>::infinity();
+
+ return func_->regularity(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 &u, VectorType const &v,
@@ -83,11 +90,12 @@ template <int dimension> class LocalNonlinearity {
// TODO: do not evaluate at zero
void addGradient(VectorType const &x, VectorType &y) const {
- if (x == VectorType(0))
+ if (x.two_norm() <
+ 1e-8 or std::isinf(func_->regularity(x.two_norm()))) // TODO
return;
VectorType tmp;
- upperGradient(x, tmp); // TODO
+ upperGradient(x, tmp); // upper and lower gradient coincide in this case
y += tmp;
}
diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index 70f70581..d571ae7b 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -83,22 +83,9 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
// determine truncation pattern
linearization.truncation.resize(u.size());
linearization.truncation.unsetAll();
- Interval<double> ab;
- for (int i = 0; i < u.size(); ++i) {
- for (int j = 0; j < block_size; ++j) {
- // TODO: handle smooth domain
- // problem.phi.smoothDomain(i, u[i][j], regularity_tol, ab, j);
-
- // // if (phi is not regular at u) or (u is not in smooth domain) or
- // (component should be ignored)
- // // truncate this component
- // if ((problem.phi.regularity(i, u[i][j], j) > regularity_tol)
- // or (u[i][j] < ab[0])
- // or (u[i][j] > ab[1])
- // or ignore[i][j])
- // linearization.truncation[i][j] = true;
- }
- }
+ for (int i = 0; i < u.size(); ++i)
+ if (problem.phi.regularity(i, u[i]) > 1e8) // TODO
+ linearization.truncation[i] = true;
// construct sparsity pattern for linearization
Dune::MatrixIndexSet indices(problem.A.N(), problem.A.M());
diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/nicefunction.hh
index 6dedd620..575b473e 100644
--- a/dune/tectonic/nicefunction.hh
+++ b/dune/tectonic/nicefunction.hh
@@ -18,6 +18,10 @@ class NiceFunction : public VirtualFunction<double, double> {
double virtual second_deriv(double x) const {
DUNE_THROW(NotImplemented, "second derivative not implemented");
}
+
+ double virtual regularity(double s) const {
+ DUNE_THROW(NotImplemented, "regularity not implemented");
+ }
};
class RuinaFunction : public NiceFunction {
@@ -82,14 +86,20 @@ class RuinaFunction : public NiceFunction {
= a/x
*/
double virtual second_deriv(double s) const {
- if (eta * s < rho)
+ // includes the case eta * s = rho for which there is no second derivative
+ if (eta * s <= rho)
return 0;
- else if (eta * s == rho)
- return 37; // TODO: not differentiable;
else
return coefficient * normalStress * (a / s);
}
+ double virtual regularity(double s) const {
+ if (eta * s == rho)
+ return std::numeric_limits<double>::infinity();
+
+ return std::abs(second_deriv(s));
+ }
+
private:
double coefficient;
double a;
@@ -116,6 +126,8 @@ class LinearFunction : public NiceFunction {
double virtual second_deriv(double) const { return 0; }
+ double virtual regularity(double s) const { return 0; }
+
private:
double coefficient;
};
@@ -153,6 +165,8 @@ class TrivialFunction : public NiceFunction {
double virtual rightDifferential(double) const { return 0; }
double virtual second_deriv(double) const { return 0; }
+
+ double virtual regularity(double) const { return 0; }
};
// slope in [n-1,n] is n
diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset
index 32366667..d5f9d7a6 100644
--- a/src/one-body-sample.parset
+++ b/src/one-body-sample.parset
@@ -49,7 +49,7 @@ normalstress = 0.1 # laursen depends a lot more on this
# http://earthquake.usgs.gov/research/physics/lab/prediction.pdf
mu = 0.5
eta = 1
-model = Ruina # TODO: Laursen doesn't work with TNNMG
+model = Ruina
[boundary.friction.ruina]
# "For rocks, typical values of A and B range from 0.005 to 0.015"
--
GitLab