ProblemStat.hpp 12.7 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
#include <dune/istl/matrix.hh>
14
15
#include <dune/grid/common/grid.hh>

16
#include <dune/amdis/AdaptInfo.hpp>
17
#include <dune/amdis/Assembler.hpp>
18
19
20
21
22
23
24
25
#include <dune/amdis/CreatorInterface.hpp>
#include <dune/amdis/DirichletBC.hpp>
#include <dune/amdis/Flag.hpp>
#include <dune/amdis/Initfile.hpp>
#include <dune/amdis/LinearAlgebra.hpp>
#include <dune/amdis/Mesh.hpp>
#include <dune/amdis/Operator.hpp>
#include <dune/amdis/ProblemStatBase.hpp>
26
#include <dune/amdis/ProblemStatTraits.hpp>
27
28
#include <dune/amdis/StandardProblemIteration.hpp>

29
#include <dune/amdis/common/OptionalDelete.hpp>
30
31
#include <dune/amdis/common/Timer.hpp>
#include <dune/amdis/common/TupleUtility.hpp>
32
#include <dune/amdis/common/TypeDefs.hpp>
33
#include <dune/amdis/common/Utility.hpp>
34

35
36
37
#include <dune/amdis/io/FileWriterInterface.hpp>

#include <dune/amdis/utility/DiscreteFunction.hpp>
38
#include <dune/amdis/utility/RangeType.hpp>
39
#include <dune/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
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
55

56
    using Element  = typename GridView::template Codim<0>::Entity;
57

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

60
    /// Dimension of the mesh
61
    static constexpr int dim = Grid::dimension;
62

63
    /// Dimension of the world
64
    static constexpr int dow = Grid::dimensionworld;
65

66
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
67

68
69
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
    using SystemVector = DOFVector<GlobalBasis, double>;
70

71
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
72

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

90
91
    /// Constructor taking additionally a reference to a mesh that is used
    /// instead of the default created mesh, \ref ProblemStat
92
    ProblemStat(std::string name, Grid* gridPtr)
93
      : ProblemStat(name)
94
    {
95
      this->grid = std::shared_ptr<Grid>(gridPtr, optional_delete(false));
96
      componentGrids.resize(nComponents, this->grid.get());
97

98
99
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
100
101
    }

102

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

113
114

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

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


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

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

139

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

146
147
148
149
150
151
152
153
    template <class Predicate, class RowTreePath, class ColTreePath>
    void addDirichletBC(Predicate const& predicate,
                        RowTreePath row, ColTreePath col,
                        double constant)
    {
      addDirichletBC(predicate, row, col, [constant](auto const&) { return constant; });
    }

154

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


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

167
168

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

Praetorius, Simon's avatar
Praetorius, Simon committed
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
203
204
205
206
207
208
209
210
211
212
  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; }


213
  public: // get-methods
214

215
    /// Returns a pointer to system-matrix, \ref systemMatrix
216
217
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
218
219


220
    /// Return a pointer to the solution, \ref solution
221
222
    std::shared_ptr<SystemVector>       getSolution()       { return solution; }
    std::shared_ptr<SystemVector const> getSolution() const { return solution; }
223

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


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


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

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

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

249
250
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
251
    void setGrid(Grid& grid_)
252
    {
253
254
      grid = std::shared_ptr<Grid>(&grid_, optional_delete(false));
      std::fill(componentGrids.begin(), componentGrids.end(), this->grid.get());
255

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


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

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

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


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

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

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

289
  protected: // initialization methods
290

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

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

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

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

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

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

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

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

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

335
    void createFileWriter();
336

Praetorius, Simon's avatar
Praetorius, Simon committed
337

338
  private:
339
340
341
342
343
344
    /// 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;
345

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

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

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

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

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

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

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

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

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


376
377
  private: // some internal data-structures

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

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


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

389
} // end namespace AMDiS
390
391
392

#include "ProblemStat.inc.hpp"