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

myblockproblem.hh

Blame
  • user avatar
    Elias Pipping authored and Elias Pipping committed
    1ee8f2a4
    History
    Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    myblockproblem.hh 4.56 KiB
    // Based on dune/tnnmg/problem-classes/blocknonlineargsproblem.hh
    
    #ifndef MY_BLOCK_PROBLEM_HH
    #define MY_BLOCK_PROBLEM_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 MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
    public:
      typedef MyConvexProblemTypeTEMPLATE MyConvexProblemType;
      typedef typename MyConvexProblemType::NonlinearityType NonlinearityType;
      typedef typename MyConvexProblemType::VectorType VectorType;
      typedef typename MyConvexProblemType::MatrixType MatrixType;
      typedef typename MyConvexProblemType::LocalVectorType LocalVectorType;
      typedef typename MyConvexProblemType::LocalMatrixType LocalMatrixType;
    
      static const int block_size = MyConvexProblemType::block_size;
    
      /** \brief Solves one local system using a scalar Gauss-Seidel method */
      class IterateObject;
    
      MyBlockProblem(MyConvexProblemType& problem) : problem(problem) {
        // TODO: Is it clever to create a bisection here?
        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
      MyConvexProblemType& problem;
    
      // commonly used minimization stuff
      Bisection bisection;
    };
    
    /** \brief Solves one local system using a scalar Gauss-Seidel method */
    template <class MyConvexProblemTypeTEMPLATE>
    class MyBlockProblem<MyConvexProblemTypeTEMPLATE>::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, MyConvexProblemType& problem)
          : problem(problem), bisection(bisection) {};
    
    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;
          Dune::minimise(localJ, ui, 10); // FIXME: hardcoded value
        }
      }
    
    private:
      // problem data
      MyConvexProblemType& problem;
    
      // commonly used minimization stuff
      Bisection bisection;
    
      // 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;
    };
    
    #endif