ProblemStat.hpp 13.1 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
  template <class Traits>
44
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
45
      : public ProblemStatBase
46
      , public StandardProblemIteration
47
  {
48
    using Self = ProblemStat;
49

50
  public: // typedefs and static constants
51

52
    using GlobalBasis = typename Traits::GlobalBasis;
53
54
55
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
    using Element  = typename GridView::template Codim<0>::Entity;
56
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
57

58
    /// Dimension of the grid
59
    static constexpr int dim = Grid::dimension;
60

61
    /// Dimension of the world
62
    static constexpr int dow = Grid::dimensionworld;
63

64
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
65
    using SystemVector = DOFVector<Traits, double>;
66

67
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
68

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

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

97

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

108
109

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

114
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
115
116
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath = {}, ColTreePath = {});
117
    /** @} */
118
119
120


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

125
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
126
127
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath = {});
128
    /** @} */
129

130

131
    /// Adds a Dirichlet boundary condition
132
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
133
    void addDirichletBC(Predicate const& predicate,
134
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
135
                        Values const& values);
136

137
  public:
138

139
    /// Implementation of \ref ProblemStatBase::solve
140
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
141
142
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
143

144
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
145
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
146
147
148
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
149

150
    /// Writes output files.
151
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
152

Praetorius, Simon's avatar
Praetorius, Simon committed
153

154
  public: // get-methods
155

156
157
158
    /// Returns a pointer to system-matrix, \ref systemMatrix_
    std::shared_ptr<SystemMatrix> getSystemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
159

160
    /// Returns a pointer to the solution vector, \ref solution_
161
    std::shared_ptr<SystemVector> getSolutionVector() const
162
    {
163
      return solution_;
164
165
166
167
168
169
170
    }

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

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


183
    /// Return a point to the rhs system-vector, \ref rhs
184
185
    std::shared_ptr<SystemVector> getRhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
186
187


188
    /// Return a pointer to the linear solver, \ref linearSolver
189
    std::shared_ptr<LinearSolverType> getSolver() { return linearSolver_; }
190

Praetorius, Simon's avatar
Praetorius, Simon committed
191
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
192
    {
193
      linearSolver_ = solver;
194
195
    }

196
    /// Return a pointer to the grid, \ref grid
197
    std::shared_ptr<Grid> getGrid() { return grid_; }
198

Praetorius, Simon's avatar
Praetorius, Simon committed
199
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
200
    /// matrices and vectors, as well as the file-writer.
201
    void setGrid(Grid& grid)
202
    {
203
204
      grid_ = Dune::stackobject_to_shared_ptr(grid);

205
      createGlobalBasis();
206
207
208
209
      createMatricesAndVectors();
      createFileWriter();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
210
211
    /// Return the gridView of the leaf-level
    auto const& gridView() { return globalBasis_->gridView(); }
212

213
    /// Return the \ref feSpaces
214
    std::shared_ptr<GlobalBasis> const& globalBasis() { return globalBasis_; }
215
216


217
    /// Implementation of \ref ProblemStatBase::getName
218
219
    virtual std::string getName() const override
    {
220
      return name_;
221
    }
222

223
  protected: // initialization methods
224

225
226
227
228
229
230
231
232
233
234
235
    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)
    {
236
      assert( bool(grid_) );
237
      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
238
      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
239
240
241
242
243
244
245
    }

    void createGlobalBasis(std::false_type)
    {
      error_exit("Can not create GlobalBasis from type. Pass a BasisCreator instead!");
    }

246
    void createGrid()
247
    {
248
249
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
250
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
251

252
      grid_ = MeshCreator<Grid>::create(gridName_);
253

254
      msg("Create grid:");
255
256
257
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
258
      msg("");
259
    }
260

261
262
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
263
264
265
266
      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{});
267
    }
268

269
270
    void createMatricesAndVectors()
    {
271
272
273
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
274

275
      AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
276
277
      {
        std::string i = to_string(treePath);
278
279
280
        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
281
282
283
      });

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

286
    void createSolver()
287
    {
288
      std::string solverName = "cg";
289
      Parameters::get(name_ + "->solver->name", solverName);
290

291
      auto solverCreator
292
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
293

294
      linearSolver_ = solverCreator->create(name_ + "->solver");
295
    }
296

297
    void createMarker();
298

299
    void createFileWriter();
300

Praetorius, Simon's avatar
Praetorius, Simon committed
301

302
303
304
305
306
307
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
  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.
334
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
335
336


337
  private:
338
    /// Name of this problem.
339
    std::string name_;
340

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

344
345
346
    /// Number of grids
    int nGrids = 1;

347
    /// Name of the grid
348
    std::string gridName_ = "none";
349

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

354
    /// A FileWriter object
355
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
356

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

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

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

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

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

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


378
379
  private: // some internal data-structures

380
381
382
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
383
  };
384

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

400

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

405
} // end namespace AMDiS
406
407
408

#include "ProblemStat.inc.hpp"