Skip to content
Snippets Groups Projects
Select Git revision
  • 1c2d7d26fec521576ae77e39bb536b659bf0087c
  • 2016-PippingKornhuberRosenauOncken default
  • 2014-Dissertation-Pipping
  • 2013-PippingSanderKornhuber
4 results

myblockproblem.hh

Blame
  • user avatar
    Elias Pipping authored and Elias Pipping committed
    1c2d7d26
    History
    Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    myblockproblem.hh 4.59 KiB
    // 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 m 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 m,
          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?)
        {
          LocalMatrixType const* localA = NULL;
          LocalVectorType localb(problem.f[m]);
    
          typename MatrixType::row_type::ConstIterator it;
          typename MatrixType::row_type::ConstIterator end = problem.A[m].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;
          }
        }
      }
    
    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;
    };