diff --git a/dune/tectonic/globalnonlinearity.hh b/dune/tectonic/globalnonlinearity.hh index 5bf0d30f7d2695db0124f9e895c3b913b83a6e07..4ea304e0e0d5a2e6b8337413cc70e44a4e9df70e 100644 --- a/dune/tectonic/globalnonlinearity.hh +++ b/dune/tectonic/globalnonlinearity.hh @@ -10,10 +10,14 @@ #include "localnonlinearity.hh" namespace Dune { -template <int dim, class VectorType = BlockVector<FieldVector<double, dim>>, - class MatrixType = BCRSMatrix<FieldMatrix<double, dim, dim>>> +template <int dim, + class VectorTypeTEMPLATE = BlockVector<FieldVector<double, dim>>, + class MatrixTypeTEMPLATE = BCRSMatrix<FieldMatrix<double, dim, dim>>> class GlobalNonlinearity { public: + typedef VectorTypeTEMPLATE VectorType; + typedef MatrixTypeTEMPLATE MatrixType; + /* Return a restriction of the outer function to the i'th node. */ @@ -27,6 +31,25 @@ class GlobalNonlinearity { } } + void directionalDomain(VectorType const &, VectorType const &, + Interval<double> &dom) const { + dom[0] = -std::numeric_limits<double>::max(); + dom[1] = std::numeric_limits<double>::max(); + } + + void directionalSubDiff(VectorType const &u, VectorType const &v, + Interval<double> &subdifferential) const { + subdifferential[0] = subdifferential[1] = 0; + for (size_t i = 0; i < u.size(); ++i) { + Interval<double> D; + auto res = restriction(i); + // TODO: is the coordinate system the right one here? + res->directionalSubDiff(u[i], v[i], D); + subdifferential[0] += D[0]; + subdifferential[1] += D[1]; + } + } + void addHessianIndices(Dune::MatrixIndexSet &indices) const { for (size_t i = 0; i < indices.rows(); ++i) indices.add(i, i); diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh index 83fbfbc22dbbdc3fead55199a330f7aec7640f7b..563437633d03e1871d158b932d67f66f15986a48 100644 --- a/dune/tectonic/myblockproblem.hh +++ b/dune/tectonic/myblockproblem.hh @@ -10,6 +10,7 @@ #include <dune/tnnmg/problem-classes/onedconvexfunction.hh> #include "localnonlinearity.hh" +#include "mydirectionalconvexfunction.hh" #include "samplefunctional.hh" /** \brief Base class for problems where each block can be solved with a @@ -71,9 +72,29 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem { double computeDampingParameter(VectorType const &u, VectorType const &projected_v) const { - // TODO: implement - return 1.0; - }; + VectorType const v = projected_v; + + VectorType tmp = problem.f; + problem.A.mmv(u, tmp); + double const localb = tmp * v; + + problem.A.mv(v, tmp); + double const localA = tmp * v; + + /* + 1/2 <A(u + hv),u + hv> - <b, u + hv> + = 1/2 <Av,v> h^2 - <b - Au, v> h + const. + + localA = <Av,v> + localb = <b - Au, v> + */ + MyDirectionalConvexFunction<Dune::GlobalNonlinearity<block_size>> psi( + localA, localb, problem.phi, u, v); + + int bisectionsteps = 0; + Bisection bisection(0.0, 1.0, 1e-12, true, 0); // TODO + return bisection.minimize(psi, 0.0, 0.0, bisectionsteps); // TODO + } void assembleTruncate(VectorType const &u, Linearization &linearization, Dune::BitSetVector<block_size> const &ignore) const {