From b0ad2aafa7e8006890a7e9ff8e025003fd014831 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Fri, 28 Oct 2011 17:41:21 +0200
Subject: [PATCH] Preliminary local problem solver (BROKEN!)
---
src/myblockproblem.hh | 45 +++++++++++++++++++++++++++++++++++++++++++
1 file changed, 45 insertions(+)
diff --git a/src/myblockproblem.hh b/src/myblockproblem.hh
index a7c535f3..1620e376 100644
--- a/src/myblockproblem.hh
+++ b/src/myblockproblem.hh
@@ -8,6 +8,10 @@
#include <dune/tnnmg/problem-classes/nonlinearity.hh>
#include <dune/tnnmg/problem-classes/onedconvexfunction.hh>
+#include "mynonlinearity.hh"
+#include "nicefunction.hh"
+#include "samplefunctional.hh"
+
/** \brief Base class for problems where each block can be solved with a scalar
* Gauss-Seidel method */
template <class ConvexProblemTypeTEMPLATE> class MyBlockProblem {
@@ -90,6 +94,47 @@ class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject {
void solveLocalProblem(
LocalVectorType& ui, int i,
const typename Dune::BitSetVector<block_size>::const_reference ignore) {
+ // Note: problem.Am and problem.a will be ignored
+ {
+ int const m = i;
+
+ LocalMatrixType const* localA = NULL;
+ LocalVectorType localb(0);
+
+ typename MatrixType::row_type::ConstIterator it;
+ typename MatrixType::row_type::ConstIterator end = problem.A[i].end();
+ for (it = problem.A[m].begin(); it != end; ++it) {
+ int j = it.index();
+ if (j == m)
+ localA = &(*it); // localA = &A[m][m]
+
+ (*it).umv(u[j], localb); // localb += A[m][j] * u[j]
+ LocalMatrixType foo = problem.A[j][m];
+ foo -= problem.A[m][j];
+ // if (foo.infinity_norm() != 0)
+ // std::cout << "Indices: " << col << ", " << m << std::endl;
+ }
+ localb -= problem.f[m];
+ assert(localA != NULL);
+
+ // FIXME: Hardcoding a fixed function here for now
+ Dune::TrivialFunction func;
+ Dune::MyNonlinearity<block_size> phi(func);
+ Dune::SampleFunctional<block_size> localJ(*localA, localb, phi);
+
+ LocalVectorType ui_copy = ui;
+
+ // std::cout << "Norm before: " << ui.two_norm() << std::endl;
+ LocalVectorType correction;
+ for (size_t i = 1; i <= 5; ++i) {
+ Dune::minimise(localJ, ui_copy, correction);
+ ui_copy += correction; // This should be correct
+ }
+
+ ui = ui_copy;
+ return;
+ }
+
const LocalMatrixType* Aii;
LocalVectorType b(0.0);
--
GitLab