From 3793cea110b7c68f91833688b3053387050388b9 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 8 May 2012 14:52:46 +0200
Subject: [PATCH] Do not use bisection for the Ruina update
---
src/compute_state_ruina.cc | 57 ++++----------------------------------
1 file changed, 5 insertions(+), 52 deletions(-)
diff --git a/src/compute_state_ruina.cc b/src/compute_state_ruina.cc
index ddba439c..a5ae9bbf 100644
--- a/src/compute_state_ruina.cc
+++ b/src/compute_state_ruina.cc
@@ -1,61 +1,14 @@
#include <cassert>
-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/fufem/interval.hh>
-
-#include <dune/tnnmg/problem-classes/bisection.hh>
-
-#include <dune/tectonic/mydirectionalconvexfunction.hh>
+#include <cmath>
#include "compute_state_ruina.hh"
-// Ruina's nonlinearity
-class Ruina {
-public:
- typedef Dune::FieldVector<double, 1> VectorType;
- typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
-
- Ruina(double h, double uol) : h(h), uol(uol) {}
-
- void directionalSubDiff(VectorType const &u, VectorType const &v,
- Interval<double> &D) const {
- if (uol == 0)
- D[0] = D[1] = v[0] * uol * (u[0]);
- else
- D[0] = D[1] = v[0] * uol * (u[0] + std::log(uol / h));
- }
-
- void directionalDomain(VectorType const &, VectorType const &,
- Interval<double> &dom) const {
- dom[0] = -std::numeric_limits<double>::max();
- dom[1] = std::numeric_limits<double>::max();
- }
-
-private:
- double const h;
- double const uol;
-};
-
-double state_update_ruina_bisection(double h, double uol, double old_state) {
- Ruina::VectorType const start(0);
- Ruina::VectorType const direction(1);
- Ruina const phi(h, uol);
- MyDirectionalConvexFunction<Ruina> const J(1.0, old_state, phi, start,
- direction);
-
- int bisectionsteps = 0;
- Bisection const bisection(0.0, 1.0, 1e-12, true, 0); // TODO
- return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO
-}
-
double state_update_ruina(double h, double uol, // unorm/L
double old_state) {
- double const uol_log_uol = (uol == 0) ? 0.0 : uol * std::log(uol);
-
- double ret1 = state_update_ruina_bisection(h, uol, old_state);
- assert(std::abs(ret1 + uol * ret1 - old_state + uol_log_uol -
- uol * std::log(h)) < 1e-14);
+ if (uol == 0)
+ return old_state;
- return ret1;
+ double const ret = old_state - uol * std::log(uol / h);
+ return ret / (1 + uol);
}
--
GitLab