Preliminary local problem solver (BROKEN!)

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);