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
    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
    /// Return a pointer to the solution, \ref solution
211
212
    std::shared_ptr<SystemVector>       getSolution()       { return solution; }
    std::shared_ptr<SystemVector const> getSolution() const { return solution; }
213

214
    /// Return a view to a solution component (as shared_ptr)
215
    template <class TreePath>
216
    auto getSolution(TreePath path) const
217
    {
218
      auto tp = makeTreePath(path);
219
      return makeDOFVectorView(*solution, tp);
220
    }
221
222


223
    /// Return a point to the rhs system-vector, \ref rhs
224
225
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
226
227


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

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

236
237
    /// Return a pointer to the grid, \ref grid
    std::shared_ptr<Grid> getGrid() { return grid; }
238

239
240
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
241
    void setGrid(Grid& grid_)
242
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
243
      grid = Dune::stackobject_to_shared_ptr(grid_);
244
      std::fill(componentGrids.begin(), componentGrids.end(), this->grid.get());
245

246
      createGlobalBasis();
247
248
249
250
251
      createMatricesAndVectors();
      createFileWriter();
    }


252
    /// Return the gridView of the leaf-level
253
    auto leafGridView() { return grid->leafGridView(); }
254
255

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

258
    /// Return the \ref feSpaces
259
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis; }
260
261


262
    /// Implementation of \ref ProblemStatBase::getName
263
264
265
266
    virtual std::string getName() const override
    {
      return name;
    }
267

268
    /// Return a vector of names of the problem components
269
270
271
272
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
273

274
    std::size_t getNumComponents() const
275
276
277
278
    {
      return nComponents;
    }

279
  protected: // initialization methods
280

281
    void createGrid()
282
    {
283
284
285
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
      test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
286

287
      grid = MeshCreator<Grid>::create(gridName);
288

289
290
291
292
      msg("Create grid:");
      msg("#elements = "   , grid->size(0));
      msg("#faces/edges = ", grid->size(1));
      msg("#vertices = "   , grid->size(dim));
293
      msg("");
294
    }
295

296
    void createGlobalBasis()
297
    {
298
      globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(leafGridView()));
299
300
301
302
      globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis->localView().tree());
      matrixOperators.init(*globalTree);
      rhsOperators.init(*globalTree);
      constraints.init(*globalTree);
303
    }
304

305
306
    void createMatricesAndVectors()
    {
307
308
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
      solution = std::make_shared<SystemVector>(*globalBasis, componentNames[0]);
309

310
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
311
      rhs = std::make_shared<SystemVector>(*globalBasis, rhsNames[0]);
312
    }
313

314
    void createSolver()
315
    {
316
317
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
318

319
320
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
321

322
      linearSolver = solverCreator->create(name + "->solver");
323
    }
324

325
    void createFileWriter();
326

Praetorius, Simon's avatar
Praetorius, Simon committed
327

328
  private:
329
330
331
332
333
334
    /// 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;
335

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

339
    /// Name of the mesh
340
    std::string gridName = "none";
341

342
    /// Pointer to the meshes for the different problem components
343
    std::vector<Grid*> componentGrids;
344

345
    /// FE spaces of this problem.
346
    std::shared_ptr<GlobalBasis> globalBasis; // eventuell const
347
    std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
348

349
    /// A FileWriter object
350
    std::list<std::shared_ptr<FileWriterInterface>> filewriter;
351

352
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
353
    std::shared_ptr<LinearSolverType> linearSolver;
354

355
    /// A block-matrix that is filled during assembling
356
    std::shared_ptr<SystemMatrix> systemMatrix;
357

358
    /// A block-vector with the solution components
359
    std::shared_ptr<SystemVector> solution;
360
361

    /// A block-vector (load-vector) corresponding to the right.hand side
362
    /// of the equation, filled during assembling
363
    std::shared_ptr<SystemVector> rhs;
364
365


366
367
  private: // some internal data-structures

368
369
    MatrixOperators<GlobalBasis> matrixOperators;
    VectorOperators<GlobalBasis> rhsOperators;
370

371
    Constraints<GlobalBasis> constraints;
372
  };
373
374


375
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
376
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
377
#endif
378

379
} // end namespace AMDiS
380
381
382

#include "ProblemStat.inc.hpp"