ProblemStat.hpp 12.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
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#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>
#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>

#include <amdis/utility/TreePath.hpp>
37

38
namespace AMDiS
39
{
40
  template <class Traits>
41
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
42
      : public ProblemStatBase
43
      , public StandardProblemIteration
44
  {
45
    using Self = ProblemStat;
46

47
  public: // typedefs and static constants
48

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

    /// Dimension of the mesh
56
    static constexpr int dim = Grid::dimension;
57

58
    /// Dimension of the world
59
    static constexpr int dow = Grid::dimensionworld;
60

61
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
62
    using SystemVector = DOFVector<Traits, double>;
63

64
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
65

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

76
77
    /// Constructor taking additionally a reference to a mesh that is used
    /// instead of the default created mesh, \ref ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
78
    ProblemStat(std::string name, Grid& grid)
79
      : ProblemStat(std::move(name))
80
    {
81
82
83
84
85
86
87
88
89
90
91
      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_);
92
93
    }

94

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

105
106

    /// Adds an operator to \ref A.
107
    /** @{ */
108
109
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath, ColTreePath,
110
                           double* factor = nullptr, double* estFactor = nullptr);
111

112
    // operator evaluated on the boundary
113
114
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath, ColTreePath,
115
                           double* factor = nullptr, double* estFactor = nullptr);
116
    /** @} */
117
118
119


    /// Adds an operator to \ref rhs.
120
    /** @{ */
121
122
    template <class Operator, class TreePath>
    void addVectorOperator(Operator const& op, TreePath,
123
                           double* factor = nullptr, double* estFactor = nullptr);
124

125
    // operator evaluated on the boundary
126
127
    template <class Operator, class TreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath,
128
                           double* factor = nullptr, double* estFactor = nullptr);
129
    /** @} */
130

131

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

138
  public:
139

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


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

152
153

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

Praetorius, Simon's avatar
Praetorius, Simon committed
156

157
  public: // get-methods
158

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

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

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

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


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


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

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

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

202
203
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
204
    void setGrid(Grid& grid)
205
    {
206
207
      grid_ = Dune::stackobject_to_shared_ptr(grid);

208
      createGlobalBasis();
209
210
211
212
213
      createMatricesAndVectors();
      createFileWriter();
    }


214
    /// Return the gridView of the leaf-level
215
216
217
218
    auto leafGridView() { return grid_->leafGridView(); }

    /// Return the gridView of levle `level`
    auto levelGridView(int level) { return grid_->levelGridView(level); }
219

220
    /// Return the \ref feSpaces
221
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { 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
    void createGrid()
233
    {
234
235
236
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
      test_exit(!gridName_.empty(), "No mesh name specified for '", name_, "->mesh'!");
237

238
      grid_ = MeshCreator<Grid>::create(gridName_);
239

240
      msg("Create grid:");
241
242
243
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
244
      msg("");
245
    }
246

247
    void createGlobalBasis()
248
    {
249
250
      globalBasis_ = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid_->leafGridView()));
      initGlobalBasis(*globalBasis_);
251
252
253
254
    }

    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
255
256
257
258
      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{});
259
    }
260

261
262
    void createMatricesAndVectors()
    {
263
264
265
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
266
    }
267

268
    void createSolver()
269
    {
270
      std::string solverName = "cg";
271
      Parameters::get(name_ + "->solver->name", solverName);
272

273
      auto solverCreator
274
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
275

276
      linearSolver_ = solverCreator->create(name_ + "->solver");
277
    }
278

279
    void createFileWriter();
280

Praetorius, Simon's avatar
Praetorius, Simon committed
281

282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
  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.
    virtual Flag markElements(AdaptInfo& adaptInfo) override { return 0; }


317
  private:
318
    /// Name of this problem.
319
    std::string name_;
320

321
    /// Grid of this problem.
322
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
323

324
    /// Name of the mesh
325
    std::string gridName_ = "none";
326

327
    /// FE spaces of this problem.
328
329
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
330

331
    /// A FileWriter object
332
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
333

334
    /// An object of the linearSolver Interface
335
    std::shared_ptr<LinearSolverType> linearSolver_;
336

337
    /// A block-matrix that is filled during assembling
338
    std::shared_ptr<SystemMatrix> systemMatrix_;
339

340
    /// A block-vector with the solution components
341
    std::shared_ptr<SystemVector> solution_;
342
343

    /// A block-vector (load-vector) corresponding to the right.hand side
344
    /// of the equation, filled during assembling
345
    std::shared_ptr<SystemVector> rhs_;
346
347


348
349
  private: // some internal data-structures

350
351
352
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
353
  };
354

355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
#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};
  }

370

371
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
372
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
373
#endif
374

375
} // end namespace AMDiS
376
377
378

#include "ProblemStat.inc.hpp"