ProblemStat.hpp 12.9 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
#include <amdis/Mesh.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>

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

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

#include <amdis/io/FileWriterInterface.hpp>

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

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

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

53
54
    friend class ProblemInstat<Traits>;

55
  public: // typedefs and static constants
56

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

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

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

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

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

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

84
85
    /// 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
86
    ProblemStat(std::string name, Grid& grid)
87
      : ProblemStat(std::move(name))
88
    {
89
90
91
92
93
94
95
96
97
98
99
      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_);
100
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
101

102

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

113
114

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

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


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

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

135

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

142
  public:
143

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

149
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
150
151
152
153
    virtual void buildAfterAdapt(AdaptInfo& adaptInfo,
                                 Flag flag,
                                 bool asmMatrix = true,
                                 bool asmVector = true) override;
154

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

Praetorius, Simon's avatar
Praetorius, Simon committed
158

159
  public: // get-methods
160

161
162
163
164
165
166
167
168
169
170
171
    /// 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_; }
172

173
174
175
176
177
178

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

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


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

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

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

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

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

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

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

    /// Return the \ref globalBasis_
    GlobalBasis& globalBasis() { return *globalBasis_; }
224
225


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

232
  protected: // initialization methods
233

234
235
236
237
238
239
240
241
242
243
244
    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)
    {
245
      assert( bool(grid_) );
246
      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
247
      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
248
249
250
251
    }

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

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

261
      grid_ = MeshCreator<Grid>::create(gridName_);
262

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

270
271
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
272
273
274
275
      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{});
276
    }
277

278
279
    void createMatricesAndVectors()
    {
280
281
282
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_);
      solution_ = std::make_shared<SystemVector>(*globalBasis_);
      rhs_ = std::make_shared<SystemVector>(*globalBasis_);
283

284
      AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
285
286
      {
        std::string i = to_string(treePath);
287
288
289
        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
290
      });
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
  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::estimate.
    virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::refineMesh.
329
    virtual Flag adaptGrid(AdaptInfo& adaptInfo) override;
330
331

    /// Implementation of \ref ProblemStatBase::markElements.
332
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
333
334


335
  private:
336
    /// Name of this problem.
337
    std::string name_;
338

339
    /// Grid of this problem.
340
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
341

342
343
344
    /// Number of grids
    int nGrids = 1;

345
    /// Name of the grid
346
    std::string gridName_ = "none";
347

348
    /// FE spaces of this problem.
349
350
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
351

352
    /// A FileWriter object
353
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
354

355
    /// Pointer to the adaptation markers
356
    std::list<std::shared_ptr<Marker<Traits> > > marker_;
357

358
    /// An object of the linearSolver Interface
359
    std::shared_ptr<LinearSolverType> linearSolver_;
360

361
    /// A block-matrix that is filled during assembling
362
    std::shared_ptr<SystemMatrix> systemMatrix_;
363

364
    /// A block-vector with the solution components
365
    std::shared_ptr<SystemVector> solution_;
366

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

371
    /// A block-vector (load-vector) corresponding to the right.hand side
372
    /// of the equation, filled during assembling
373
    std::shared_ptr<SystemVector> rhs_;
374
375


376
377
  private: // some internal data-structures

378
379
380
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
381
  };
382

383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
#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};
  }

398

399
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
400
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
401
#endif
402

403
} // end namespace AMDiS
404
405
406

#include "ProblemStat.inc.hpp"