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

53
    using Element  = typename GridView::template Codim<0>::Entity;
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 WorldVector = typename Element::Geometry::GlobalCoordinate;
62

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

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

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

78
79
    /// 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
80
    ProblemStat(std::string name, Grid& grid)
81
      : ProblemStat(name)
82
    {
83
84
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
85
86
    }

87

88
    /**
89
     * \brief Initialisation of the problem.
90
     *
91
92
93
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
94
    void initialize(Flag initFlag,
95
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
96
                    Flag adoptFlag = INIT_NOTHING);
97

98
99

    /// Adds an operator to \ref A.
100
    /** @{ */
101
102
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath, ColTreePath,
103
                           double* factor = nullptr, double* estFactor = nullptr);
104

105
    // operator evaluated on the boundary
106
107
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath, ColTreePath,
108
                           double* factor = nullptr, double* estFactor = nullptr);
109
    /** @} */
110
111
112


    /// Adds an operator to \ref rhs.
113
    /** @{ */
114
115
    template <class Operator, class TreePath>
    void addVectorOperator(Operator const& op, TreePath,
116
                           double* factor = nullptr, double* estFactor = nullptr);
117

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

124

125
    /// Adds a Dirichlet boundary condition
126
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
127
    void addDirichletBC(Predicate const& predicate,
128
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
129
                        Values const& values);
130

131
  public:
132

133
    /// Implementation of \ref ProblemStatBase::solve
134
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
135
136
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
137
138


139
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
140
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
141
142
143
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
144

145
146

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

Praetorius, Simon's avatar
Praetorius, Simon committed
149

150
151
152
153
154
155
156
157
158
159
160
161
162
  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
    {
163
164
165
166
167
      // for (std::size_t i = 0; i < getNumComponents(); ++i)
      //   if (adaptInfo.spaceToleranceReached(i))
      //     adaptInfo.allowRefinement(false, i);
      //   else
      //     adaptInfo.allowRefinement(true, i);
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190

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


191
  public: // get-methods
192

193
    /// Returns a pointer to system-matrix, \ref systemMatrix
194
195
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
196
197


198
    std::shared_ptr<SystemVector> getSolutionVector() const
199
200
201
202
203
204
205
206
207
208
209
    {
      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);
    }
210

211
212
213
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
214
    {
215
      auto&& tp = makeTreePath(path);
216
      return makeDOFVectorView(*solution, tp);
217
    }
218
219


220
    /// Return a point to the rhs system-vector, \ref rhs
221
222
    auto getRhsVector()       { return rhs; }
    auto getRhsVector() const { return rhs; }
223
224


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

Praetorius, Simon's avatar
Praetorius, Simon committed
228
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
229
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
230
      linearSolver = solver;
231
232
    }

233
234
    /// Return a pointer to the grid, \ref grid
    std::shared_ptr<Grid> getGrid() { return grid; }
235

236
237
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
238
    void setGrid(Grid& grid_)
239
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
240
      grid = Dune::stackobject_to_shared_ptr(grid_);
241
      createGlobalBasis();
242
243
244
245
      createMatricesAndVectors();
      createFileWriter();
    }

246
/*
247
    /// Return the gridView of the leaf-level
248
    auto leafGridView() { return grid->leafGridView(); }
249

250
251
252
    /// Return the gridView of levle `level`
    auto levelGridView(int level) { return grid->levelGridView(level); }
*/
253
    /// Return the \ref feSpaces
254
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis; }
255
256


257
    /// Implementation of \ref ProblemStatBase::getName
258
259
260
261
    virtual std::string getName() const override
    {
      return name;
    }
262

263
  protected: // initialization methods
264

265
    void createGrid()
266
    {
267
268
269
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
      test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
270

271
      grid = MeshCreator<Grid>::create(gridName);
272

273
274
275
276
      msg("Create grid:");
      msg("#elements = "   , grid->size(0));
      msg("#faces/edges = ", grid->size(1));
      msg("#vertices = "   , grid->size(dim));
277
      msg("");
278
    }
279

280
    void createGlobalBasis()
281
    {
282
      globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid->leafGridView()));
283
284
285
286
      globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis->localView().tree());
      matrixOperators.init(*globalTree);
      rhsOperators.init(*globalTree);
      constraints.init(*globalTree);
287
    }
288

289
290
    void createMatricesAndVectors()
    {
291
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
292
293
      solution = std::make_shared<SystemVector>(*globalBasis, "solution");
      rhs = std::make_shared<SystemVector>(*globalBasis, "rhs");
294
    }
295

296
    void createSolver()
297
    {
298
299
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
300

301
302
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
303

304
      linearSolver = solverCreator->create(name + "->solver");
305
    }
306

307
    void createFileWriter();
308

Praetorius, Simon's avatar
Praetorius, Simon committed
309

310
  private:
311
312
313
    /// Name of this problem.
    std::string name;

314
315
    /// Grid of this problem.
    std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
316

317
    /// Name of the mesh
318
    std::string gridName = "none";
319

320
    /// FE spaces of this problem.
321
    std::shared_ptr<GlobalBasis> globalBasis;
322
    std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
323

324
    /// A FileWriter object
325
    std::list<std::shared_ptr<FileWriterInterface>> filewriter;
326

327
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
328
    std::shared_ptr<LinearSolverType> linearSolver;
329

330
    /// A block-matrix that is filled during assembling
331
    std::shared_ptr<SystemMatrix> systemMatrix;
332

333
    /// A block-vector with the solution components
334
    std::shared_ptr<SystemVector> solution;
335
336

    /// A block-vector (load-vector) corresponding to the right.hand side
337
    /// of the equation, filled during assembling
338
    std::shared_ptr<SystemVector> rhs;
339
340


341
342
  private: // some internal data-structures

343
344
    MatrixOperators<GlobalBasis> matrixOperators;
    VectorOperators<GlobalBasis> rhsOperators;
345

346
    Constraints<GlobalBasis> constraints;
347
  };
348
349


350
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
351
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
352
#endif
353

354
} // end namespace AMDiS
355
356
357

#include "ProblemStat.inc.hpp"