myblockproblem.hh

// Based on dune/tnnmg/problem-classes/blocknonlineargsproblem.hh

// #include <dune/common/parametertree.hh>
#include <dune/common/bitsetvector.hh>

#include <dune/tnnmg/problem-classes/bisection.hh>
#include <dune/tnnmg/problem-classes/convexproblem.hh>
#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 {
public:
  typedef ConvexProblemTypeTEMPLATE ConvexProblemType;
  typedef typename ConvexProblemType::NonlinearityType NonlinearityType;
  typedef typename ConvexProblemType::VectorType VectorType;
  typedef typename ConvexProblemType::MatrixType MatrixType;
  typedef typename ConvexProblemType::LocalVectorType LocalVectorType;
  typedef typename ConvexProblemType::LocalMatrixType LocalMatrixType;

  static const int block_size = ConvexProblemType::block_size;

  /** \brief Solves one local system using a scalar Gauss-Seidel method */
  class IterateObject;

  MyBlockProblem(ConvexProblemType& problem) : problem(problem) {
    bisection = Bisection(0.0, 1.0, 1e-15, true, 1e-14);
  };

  /** \brief Constructs and returns an iterate object */
  IterateObject getIterateObject() { return IterateObject(bisection, problem); }

private:
  // problem data
  ConvexProblemType& problem;

  // commonly used minimization stuff
  Bisection bisection;
};

/** \brief Solves one local system using a scalar Gauss-Seidel method */
template <class ConvexProblemTypeTEMPLATE>
class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject {
  friend class MyBlockProblem;

protected:
  /** \brief Constructor, protected so only friends can instantiate it
   * \param bisection The class used to do a scalar bisection
   * \param problem The problem including quadratic part and nonlinear/nonsmooth
   * part
   */
  IterateObject(const Bisection& bisection, ConvexProblemType& problem)
      : problem(problem),
        bisection(bisection),
        local_J(1.0, 0.0, problem.phi, 0, 0) {};

public:
  /** \brief Set the current iterate */
  void setIterate(VectorType& u) {
    // TODO How should the rang-1 term be handled
    // ????????????????????????????????
    // s = problem.Am*u;
    // problem.phi.setVector(u);

    this->u = u;
    return;
  };

  /** \brief Update the i-th block of the current iterate */
  void updateIterate(const LocalVectorType& ui, int i) {
    // TODO How should the rang-1 term be handled
    // ????????????????????????????????
    // s += (ui-u[i]) * problem.Am[i];

    // for(size_t j=0; j<block_size; ++j)
    //   problem.phi.updateEntry(i, ui[j], j);

    u[i] = ui;
    return;
  };

  /** \brief Minimize a local problem using a scalar Gauss-Seidel method
   * \param[out] ui The solution
   * \param i Block number
   * \param ignore Set of degrees of freedom to leave untouched
   *
   * \return The minimizer of the local functional in the variable ui
   */
  void solveLocalProblem(
      LocalVectorType& ui, int i,
      const typename Dune::BitSetVector<block_size>::const_reference ignore) {
    // Note: problem.Am and problem.a will be ignored
    // Note: ignore is currently ignored (what's it used for anyway?)
    {
      int const m = i;

      LocalMatrixType const* localA = NULL;
      LocalVectorType localb(problem.f[m]);

      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 const j = it.index();
        if (j == m)
          localA = &(*it); // localA = &A[m][m]
        else
          it->mmv(u[j], localb); // localb += A[m][j] * u[j]
      }
      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 correction;
      for (size_t i = 1; i <= 10; ++i) { // FIXME: hardcoded value
        Dune::minimise(localJ, ui, correction);
        ui += correction;
      }
      return;
    }
  }

private:
  // problem data
  ConvexProblemType& problem;

  // commonly used minimization stuff
  Bisection bisection;
  OneDConvexFunction<NonlinearityType> local_J;

  // state data for smoothing procedure used by:
  // setIterate, updateIterate, solveLocalProblem
  VectorType u;

  /** \brief Keeps track of the total number of bisection steps that were
   * performed */
  int bisectionsteps;
};