ProblemStat.hpp 13 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:
Praetorius, Simon's avatar
Praetorius, Simon committed
70
#if 0
71
72
    /**
     * \brief Constructor. Takes the name of the problem that is used to
73
     * access values corresponding to this problem in the parameter file.
74
     **/
75
    explicit ProblemStat(std::string name)
76
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
Praetorius, Simon's avatar
Praetorius, Simon committed
77
      , name(std::move(name))
78
    {}
79

80
81
    /// 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
82
    ProblemStat(std::string name, Grid& grid)
83
      : ProblemStat(std::move(name))
84
    {
85
86
87
88
89
90
91
92
93
94
95
      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_);
96
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#endif

    /// \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_);
    }
111

112

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

123
124

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

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


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

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

145

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

152
  public:
153

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
168

169
  public: // get-methods
170

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

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

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

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


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


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

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

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

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

221
      createGlobalBasis();
222
223
224
      createMatricesAndVectors();
      createFileWriter();
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
225
#endif
226

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

230
    /// Return the \ref feSpaces
231
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
232
233


234
    /// Implementation of \ref ProblemStatBase::getName
235
236
    virtual std::string getName() const override
    {
237
      return name_;
238
    }
239

240
  protected: // initialization methods
241

242
    void createGrid()
243
    {
244
245
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
246
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
247

248
      grid_ = MeshCreator<Grid>::create(gridName_);
249

250
      msg("Create grid:");
251
252
253
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
254
      msg("");
255
    }
256

257
258
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
259
      localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
Praetorius, Simon's avatar
Praetorius, Simon committed
260
261
262
      matrixOperators.init(localView_->tree(), tag::store{});
      rhsOperators.init(localView_->tree(), tag::store{});
      constraints.init(localView_->tree(), tag::store{});
263
    }
264

265
266
    void createMatricesAndVectors()
    {
267
268
269
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
270

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

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

282
    void createSolver()
283
    {
284
      std::string solverName = "cg";
285
      Parameters::get(name_ + "->solver->name", solverName);
286

287
      auto solverCreator
288
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
289

290
      linearSolver_ = solverCreator->create(name_ + "->solver");
291
    }
292

293
    void createMarker();
294

295
    void createFileWriter();
296

Praetorius, Simon's avatar
Praetorius, Simon committed
297

298
299
300
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
  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.
330
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
331
332


333
  private:
334
    /// Name of this problem.
335
    std::string name_;
336

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

340
341
342
    /// Number of grids
    int nGrids = 1;

343
    /// Name of the grid
344
    std::string gridName_ = "none";
345

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

350
    /// A FileWriter object
351
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
352

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

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

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

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

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

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


374
375
  private: // some internal data-structures

376
377
378
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
379
  };
380

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

396

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

401
} // end namespace AMDiS
402
403
404

#include "ProblemStat.inc.hpp"