ProblemStat.hpp 13.6 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))
Praetorius, Simon's avatar
Praetorius, Simon committed
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
100
101
102
103
104
105
106
107
108

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
    ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      , name(std::move(name))
    {
      gridName = "";
      Parameters::get(name + "->mesh", gridName);

      globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
      initGlobalBasis(*globalBasis_);
    }
109

110

111
    /**
112
     * \brief Initialisation of the problem.
113
     *
114
115
116
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
117
    void initialize(Flag initFlag,
118
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
119
                    Flag adoptFlag = INIT_NOTHING);
120

121
122

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

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


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

138
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
139
140
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath = {});
141
    /** @} */
142

143

144
    /// Adds a Dirichlet boundary condition
145
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
146
    void addDirichletBC(Predicate const& predicate,
147
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
148
                        Values const& values);
149

150
  public:
151

152
    /// Implementation of \ref ProblemStatBase::solve
153
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
154
155
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
156

157
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
158
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
159
160
161
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
162

163
    /// Writes output files.
164
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
165

Praetorius, Simon's avatar
Praetorius, Simon committed
166

167
  public: // get-methods
168

169
170
171
    /// Returns a pointer to system-matrix, \ref systemMatrix_
    std::shared_ptr<SystemMatrix> getSystemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
172

173
    /// Returns a pointer to the solution vector, \ref solution_
174
    std::shared_ptr<SystemVector> getSolutionVector() const
175
    {
176
      return solution_;
177
178
179
180
181
182
183
    }

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

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


196
    /// Return a point to the rhs system-vector, \ref rhs
197
198
    std::shared_ptr<SystemVector> getRhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
199
200


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

Praetorius, Simon's avatar
Praetorius, Simon committed
204
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
205
    {
206
      linearSolver_ = solver;
207
208
    }

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

Praetorius, Simon's avatar
Praetorius, Simon committed
212
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
213
    /// matrices and vectors, as well as the file-writer.
214
    void setGrid(Grid& grid)
215
    {
216
217
      grid_ = Dune::stackobject_to_shared_ptr(grid);

218
      createGlobalBasis();
219
220
221
222
      createMatricesAndVectors();
      createFileWriter();
    }

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

226
    /// Return the \ref feSpaces
227
    std::shared_ptr<GlobalBasis> const& globalBasis() { return globalBasis_; }
228
229


230
    /// Implementation of \ref ProblemStatBase::getName
231
232
    virtual std::string getName() const override
    {
233
      return name_;
234
    }
235

236
  protected: // initialization methods
237

238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
    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)
    {
      assert( bool(grid) );
      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid->leafGridView()));
    }

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

259
    void createGrid()
260
    {
261
262
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
263
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
264

265
      grid_ = MeshCreator<Grid>::create(gridName_);
266

267
      msg("Create grid:");
268
269
270
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
271
      msg("");
272
    }
273

274
275
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
276
      localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
Praetorius, Simon's avatar
Praetorius, Simon committed
277
278
279
      matrixOperators.init(localView_->tree(), tag::store{});
      rhsOperators.init(localView_->tree(), tag::store{});
      constraints.init(localView_->tree(), tag::store{});
280
    }
281

282
283
    void createMatricesAndVectors()
    {
284
285
286
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
287

288
      AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
289
290
      {
        std::string i = to_string(treePath);
291
292
293
        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
294
295
296
      });

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

299
    void createSolver()
300
    {
301
      std::string solverName = "cg";
302
      Parameters::get(name_ + "->solver->name", solverName);
303

304
      auto solverCreator
305
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
306

307
      linearSolver_ = solverCreator->create(name_ + "->solver");
308
    }
309

310
    void createMarker();
311

312
    void createFileWriter();
313

Praetorius, Simon's avatar
Praetorius, Simon committed
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
341
342
343
344
345
346
  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.
347
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
348
349


350
  private:
351
    /// Name of this problem.
352
    std::string name_;
353

354
    /// Grid of this problem.
355
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
356

357
358
359
    /// Number of grids
    int nGrids = 1;

360
    /// Name of the grid
361
    std::string gridName_ = "none";
362

363
    /// FE spaces of this problem.
364
365
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
366

367
    /// A FileWriter object
368
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
369

370
    /// Pointer to the adaptation markers
371
    std::list<std::shared_ptr<Marker<Traits> > > marker_;
372

373
    /// An object of the linearSolver Interface
374
    std::shared_ptr<LinearSolverType> linearSolver_;
375

376
    /// A block-matrix that is filled during assembling
377
    std::shared_ptr<SystemMatrix> systemMatrix_;
378

379
    /// A block-vector with the solution components
380
    std::shared_ptr<SystemVector> solution_;
381

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

386
    /// A block-vector (load-vector) corresponding to the right.hand side
387
    /// of the equation, filled during assembling
388
    std::shared_ptr<SystemVector> rhs_;
389
390


391
392
  private: // some internal data-structures

393
394
395
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
396
  };
397

398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
#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};
  }

413

414
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
415
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
416
#endif
417

418
} // end namespace AMDiS
419
420
421

#include "ProblemStat.inc.hpp"