ProblemStat.hpp 12.2 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
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
52

53
    using Element  = typename GridView::template Codim<0>::Entity;
54

55
    static const std::size_t nComponents = 1; // TODO: count leaf nodes in GlobalBasis
56

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

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

63
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
64

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

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

70
  public:
71
72
    /**
     * \brief Constructor. Takes the name of the problem that is used to
73
     * access values correpsonding to this püroblem in the parameter file.
74
     *
75
     * Parameters read by ProblemStat, with name 'PROB'
76
     *   PROB->names: \ref componentNames
77
     **/
78
    explicit ProblemStat(std::string name)
79
80
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      , name(name)
81
    {
82
      Parameters::get(name + "->names", componentNames);
83
      for (std::size_t i = componentNames.size(); i < nComponents; ++i)
84
        componentNames.push_back("solution[" + std::to_string(i) + "]");
85
    }
86

87
88
    /// 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
89
    ProblemStat(std::string name, Grid& grid)
90
      : ProblemStat(name)
91
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
92
      this->grid = Dune::stackobject_to_shared_ptr(grid);
93
      componentGrids.resize(nComponents, this->grid.get());
94

95
96
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
97
98
    }

99

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

110
111

    /// Adds an operator to \ref A.
112
    /** @{ */
113
114
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath, ColTreePath,
115
                           double* factor = nullptr, double* estFactor = nullptr);
116

117
    // operator evaluated on the boundary
118
119
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath, ColTreePath,
120
                           double* factor = nullptr, double* estFactor = nullptr);
121
    /** @} */
122
123
124


    /// Adds an operator to \ref rhs.
125
    /** @{ */
126
127
    template <class Operator, class TreePath>
    void addVectorOperator(Operator const& op, TreePath,
128
                           double* factor = nullptr, double* estFactor = nullptr);
129

130
    // operator evaluated on the boundary
131
132
    template <class Operator, class TreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath,
133
                           double* factor = nullptr, double* estFactor = nullptr);
134
    /** @} */
135

136

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

143
  public:
144

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


151
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
152
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
153
154
155
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
156

157
158

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

Praetorius, Simon's avatar
Praetorius, Simon committed
161

162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
  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
    {
      for (std::size_t i = 0; i < getNumComponents(); ++i)
        if (adaptInfo.spaceToleranceReached(i))
          adaptInfo.allowRefinement(false, i);
        else
          adaptInfo.allowRefinement(true, i);

      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; }


203
  public: // get-methods
204

205
    /// Returns a pointer to system-matrix, \ref systemMatrix
206
207
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
208
209


210
    std::shared_ptr<SystemVector> getSolutionVector() const
211
212
213
214
215
216
217
218
219
220
221
    {
      return solution;
    }

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

223
224
225
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
226
    {
227
      auto&& tp = makeTreePath(path);
228
      return makeDOFVectorView(*solution, tp);
229
    }
230
231


232
    /// Return a point to the rhs system-vector, \ref rhs
233
234
    auto getRhsVector()       { return rhs; }
    auto getRhsVector() const { return rhs; }
235
236


237
238
    /// Return a pointer to the linear solver, \ref linearSolver
    std::shared_ptr<LinearSolverType> getSolver() { return linearSolver; }
239

Praetorius, Simon's avatar
Praetorius, Simon committed
240
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
241
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
242
      linearSolver = solver;
243
244
    }

245
246
    /// Return a pointer to the grid, \ref grid
    std::shared_ptr<Grid> getGrid() { return grid; }
247

248
249
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
250
    void setGrid(Grid& grid_)
251
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
252
      grid = Dune::stackobject_to_shared_ptr(grid_);
253
      std::fill(componentGrids.begin(), componentGrids.end(), this->grid.get());
254

255
      createGlobalBasis();
256
257
258
259
260
      createMatricesAndVectors();
      createFileWriter();
    }


261
    /// Return the gridView of the leaf-level
262
    auto leafGridView() { return grid->leafGridView(); }
263
264

    /// Return the gridView of levle `level`
265
    auto levelGridView(int level) { return grid->levelGridView(level); }
266

267
    /// Return the \ref feSpaces
268
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis; }
269
270


271
    /// Implementation of \ref ProblemStatBase::getName
272
273
274
275
    virtual std::string getName() const override
    {
      return name;
    }
276

277
    /// Return a vector of names of the problem components
278
279
280
281
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
282

283
    std::size_t getNumComponents() const
284
285
286
287
    {
      return nComponents;
    }

288
  protected: // initialization methods
289

290
    void createGrid()
291
    {
292
293
294
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
      test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
295

296
      grid = MeshCreator<Grid>::create(gridName);
297

298
299
300
301
      msg("Create grid:");
      msg("#elements = "   , grid->size(0));
      msg("#faces/edges = ", grid->size(1));
      msg("#vertices = "   , grid->size(dim));
302
      msg("");
303
    }
304

305
    void createGlobalBasis()
306
    {
307
      globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(leafGridView()));
308
309
310
311
      globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis->localView().tree());
      matrixOperators.init(*globalTree);
      rhsOperators.init(*globalTree);
      constraints.init(*globalTree);
312
    }
313

314
315
    void createMatricesAndVectors()
    {
316
317
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
      solution = std::make_shared<SystemVector>(*globalBasis, componentNames[0]);
318

319
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
320
      rhs = std::make_shared<SystemVector>(*globalBasis, rhsNames[0]);
321
    }
322

323
    void createSolver()
324
    {
325
326
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
327

328
329
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
330

331
      linearSolver = solverCreator->create(name + "->solver");
332
    }
333

334
    void createFileWriter();
335

Praetorius, Simon's avatar
Praetorius, Simon committed
336

337
  private:
338
339
340
341
342
343
    /// Name of this problem.
    std::string name;

    /// Stores the names for all components. Is used for naming the solution
    /// vectors, \ref solution.
    std::vector<std::string> componentNames;
344

345
346
    /// Grid of this problem.
    std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
347

348
    /// Name of the mesh
349
    std::string gridName = "none";
350

351
    /// Pointer to the meshes for the different problem components
352
    std::vector<Grid*> componentGrids;
353

354
    /// FE spaces of this problem.
355
    std::shared_ptr<GlobalBasis> globalBasis;
356
    std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
357

358
    /// A FileWriter object
359
    std::list<std::shared_ptr<FileWriterInterface>> filewriter;
360

361
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
362
    std::shared_ptr<LinearSolverType> linearSolver;
363

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

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

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


375
376
  private: // some internal data-structures

377
378
    MatrixOperators<GlobalBasis> matrixOperators;
    VectorOperators<GlobalBasis> rhsOperators;
379

380
    Constraints<GlobalBasis> constraints;
381
  };
382
383


384
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
385
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
386
#endif
387

388
} // end namespace AMDiS
389
390
391

#include "ProblemStat.inc.hpp"