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;