ProblemStat.hpp 13.4 KB
Newer Older
1 2 3 4
#pragma once

#include <list>
#include <map>
5 6 7 8 9
#include <memory>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
10

11 12
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
13 14
#include <dune/grid/common/grid.hh>

15 16 17 18 19 20 21
#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/DirichletBC.hpp>
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
22
#include <amdis/LocalAssemblerList.hpp>
23
#include <amdis/Marker.hpp>
24 25 26 27 28 29 30 31 32 33 34 35 36 37
#include <amdis/Mesh.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>

#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/TypeDefs.hpp>
#include <amdis/common/Utility.hpp>

#include <amdis/GridFunctions.hpp>
#include <amdis/gridfunctions/DOFVectorView.hpp>

#include <amdis/io/FileWriterInterface.hpp>

38
#include <amdis/utility/TreeData.hpp>
39
#include <amdis/utility/TreePath.hpp>
40

41
namespace AMDiS
42
{
43 44 45 46
  // forward declaration
  template <class Traits>
  class ProblemInstat;

47
  template <class Traits>
48
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
49
      : public ProblemStatBase
50
      , public StandardProblemIteration
51
  {
52
    using Self = ProblemStat;
53

54 55
    friend class ProblemInstat<Traits>;

56
  public: // typedefs and static constants
57

58
    using GlobalBasis = typename Traits::GlobalBasis;
59 60 61
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
    using Element  = typename GridView::template Codim<0>::Entity;
62
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
63

64
    /// Dimension of the grid
65
    static constexpr int dim = Grid::dimension;
66

67
    /// Dimension of the world
68
    static constexpr int dow = Grid::dimensionworld;
69

70
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
71
    using SystemVector = DOFVector<Traits, double>;
72

73
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
74

75
  public:
76 77
    /**
     * \brief Constructor. Takes the name of the problem that is used to
78
     * access values corresponding to this problem in the parameter file.
79
     **/
80
    explicit ProblemStat(std::string name)
81
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
82
      , name_(std::move(name))
83
    {}
84

85 86
    /// Constructor taking additionally a reference to a grid that is used
    /// instead of the default created grid, \ref ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
87
    ProblemStat(std::string name, Grid& grid)
88
      : ProblemStat(std::move(name))
89
    {
90 91 92 93 94 95 96 97 98 99 100
      grid_ = Dune::stackobject_to_shared_ptr(grid);
      Parameters::get(name_ + "->mesh", gridName_);
    }

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
    ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
      : ProblemStat(std::move(name), grid)
    {
      globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
      initGlobalBasis(*globalBasis_);
101
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
102

103

104
    /**
105
     * \brief Initialisation of the problem.
106
     *
107 108 109
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
110
    void initialize(Flag initFlag,
111
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
112
                    Flag adoptFlag = INIT_NOTHING);
113

114 115

    /// Adds an operator to \ref A.
116
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
117 118
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath = {}, ColTreePath = {});
119

120
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
121 122
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath = {}, ColTreePath = {});
123
    /** @} */
124 125 126


    /// Adds an operator to \ref rhs.
127
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
128 129
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(Operator const& op, TreePath = {});
130

131
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
132 133
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath = {});
134
    /** @} */
135

136

137
    /// Adds a Dirichlet boundary condition
138
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
139
    void addDirichletBC(Predicate const& predicate,
140
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
141
                        Values const& values);
142

143
  public:
144

145
    /// Implementation of \ref ProblemStatBase::solve
146
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
147 148
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
149

150
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
151
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
152 153 154
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
155

156
    /// Writes output files.
157
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
158

Praetorius, Simon's avatar
Praetorius, Simon committed
159

160
  public: // get-methods
161

162 163 164 165 166 167 168 169 170 171 172
    /// Returns a reference to system-matrix, \ref systemMatrix_
    SystemMatrix&       getSystemMatrix()       { return *systemMatrix_; }
    SystemMatrix const& getSystemMatrix() const { return *systemMatrix_; }

    /// Returns a reference to the solution vector, \ref solution_
    SystemVector&       getSolutionVector()       { return *solution_; }
    SystemVector const& getSolutionVector() const { return *solution_; }

    /// Return a reference to the rhs system-vector, \ref rhs
    SystemVector&       getRhsVector()       { return *rhs_; }
    SystemVector const& getRhsVector() const { return *rhs_; }
173

174 175 176 177 178 179

    /// Return a mutable view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {})
    {
      auto&& tp = makeTreePath(path);
180
      return makeDOFVectorView(*solution_, tp);
181
    }
182

183 184 185
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
186
    {
187
      auto&& tp = makeTreePath(path);
188
      return makeDOFVectorView(*solution_, tp);
189
    }
190 191


192 193 194
    /// Return a reference to the linear solver, \ref linearSolver
    LinearSolverType&       getSolver()       { return *linearSolver_; }
    LinearSolverType const& getSolver() const { return *linearSolver_; }
195

196
    /// Set a new linear solver for the problem
Praetorius, Simon's avatar
Praetorius, Simon committed
197
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
198
    {
199
      linearSolver_ = solver;
200 201
    }

202 203 204
    /// Return a reference to the grid, \ref grid
    Grid&       grid()       { return *grid_; }
    Grid const& grid() const { return *grid_; }
205

206
    /// Set the grid. Stores pointer and initializes feSpaces
207
    /// matrices and vectors, as well as the file-writer.
208
    void setGrid(std::shared_ptr<Grid> const& grid)
209
    {
210
      grid_ = grid;
211

212
      createGlobalBasis();
213 214 215 216
      createMatricesAndVectors();
      createFileWriter();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
217
    /// Return the gridView of the leaf-level
218
    GridView const& gridView() { return globalBasis_->gridView(); }
219

220
    /// Return the \ref feSpaces
221
    GlobalBasis const& globalBasis() { return *globalBasis_; }
222 223


224
    /// Implementation of \ref ProblemStatBase::getName
225 226
    virtual std::string getName() const override
    {
227
      return name_;
228
    }
229

230
  protected: // initialization methods
231

232 233 234 235 236 237 238 239 240 241 242
    template <class T, class GV>
    using HasCreate = decltype(T::create(std::declval<GV>()));

    void createGlobalBasis()
    {
      createGlobalBasis(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
      initGlobalBasis(*globalBasis_);
    }

    void createGlobalBasis(std::true_type)
    {
243
      assert( bool(grid_) );
244
      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
245
      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
246 247 248 249
    }

    void createGlobalBasis(std::false_type)
    {
250
      error_exit("Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
251 252
    }

253
    void createGrid()
254
    {
255 256
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
257
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
258

259
      grid_ = MeshCreator<Grid>::create(gridName_);
260

261
      msg("Create grid:");
262 263 264
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
265
      msg("");
266
    }
267

268 269
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
270 271 272 273
      localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis_->localView());
      matrixOperators_.init(localView_->tree(), tag::store{});
      rhsOperators_.init(localView_->tree(), tag::store{});
      constraints_.init(localView_->tree(), tag::store{});
274
    }
275

276 277
    void createMatricesAndVectors()
    {
278 279 280
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
281

282
      AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
283 284
      {
        std::string i = to_string(treePath);
285 286 287
        estimates_[i].resize(globalBasis_->gridView().indexSet().size(0));
        for (std::size_t j = 0; j < estimates_[i].size(); j++)
          estimates_[i][j] = 0.0; // TODO: Remove when estimate() is implemented
288 289 290
      });

      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
291
    }
292

293
    void createSolver()
294
    {
295
      std::string solverName = "cg";
296
      Parameters::get(name_ + "->solver->name", solverName);
297

298
      auto solverCreator
299
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
300

301
      linearSolver_ = solverCreator->create(name_ + "->solver");
302
    }
303

304
    void createMarker();
305

306
    void createFileWriter();
307

Praetorius, Simon's avatar
Praetorius, Simon committed
308

309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340
  public: // implementation of iteration interface methods

    /**
      * \brief Determines the execution order of the single adaption steps.
      *
      * If adapt is true, mesh adaption will be performed. This allows to avoid
      * mesh adaption, e.g. in timestep adaption loops of timestep adaptive
      * strategies.
      *
      * Implementation of \ref StandardProblemIteration::oneIteration.
      **/
    virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
    {
      return StandardProblemIteration::oneIteration(adaptInfo, toDo);
    }

    /// Implementation of \ref ProblemStatBase::buildBeforeRefine.
    virtual void buildBeforeRefine(AdaptInfo&, Flag) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::buildBeforeCoarsen.
    virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::estimate.
    virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::refineMesh.
    virtual Flag refineMesh(AdaptInfo& adaptInfo) override { return 0; }

    /// Implementation of \ref ProblemStatBase::coarsenMesh.
    virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return 0; }

    /// Implementation of \ref ProblemStatBase::markElements.
341
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
342 343


344
  private:
345
    /// Name of this problem.
346
    std::string name_;
347

348
    /// Grid of this problem.
349
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
350

351 352 353
    /// Number of grids
    int nGrids = 1;

354
    /// Name of the grid
355
    std::string gridName_ = "none";
356

357
    /// FE spaces of this problem.
358 359
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
360

361
    /// A FileWriter object
362
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
363

364
    /// Pointer to the adaptation markers
365
    std::list<std::shared_ptr<Marker<Traits> > > marker_;
366

367
    /// An object of the linearSolver Interface
368
    std::shared_ptr<LinearSolverType> linearSolver_;
369

370
    /// A block-matrix that is filled during assembling
371
    std::shared_ptr<SystemMatrix> systemMatrix_;
372

373
    /// A block-vector with the solution components
374
    std::shared_ptr<SystemVector> solution_;
375

376
    /// A vector with the local element error estimates
377 378
    /// for each node in the basis tree, indexed by [to_string(treePath)][element index]
    std::map<std::string, std::vector<double> > estimates_;
379

380
    /// A block-vector (load-vector) corresponding to the right.hand side
381
    /// of the equation, filled during assembling
382
    std::shared_ptr<SystemVector> rhs_;
383 384


385 386
  private: // some internal data-structures

387 388 389
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
390
  };
391

392 393 394 395 396 397 398 399 400 401 402 403 404 405 406
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
  // Deduction rule
  template <class Grid, class GlobalBasis>
  ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
    -> ProblemStat<DefaultProblemTraits<GlobalBasis>>;
#endif

  // Generator for ProblemStat with given Grid and GlobalBasis
  template <class Grid, class GlobalBasis>
  ProblemStat<DefaultProblemTraits<GlobalBasis>>
  makeProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
  {
    return {std::move(name), grid, globalBasis};
  }

407

408
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
409
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
410
#endif
411

412
} // end namespace AMDiS
413 414 415

#include "ProblemStat.inc.hpp"