ProblemStat.hpp 13.6 KB
Newer Older
1
2
#pragma once

3
4
5
6
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

7
8
9
10
11
12
#include <tuple>
#include <string>
#include <memory>
#include <list>
#include <map>

13
14
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
15
#include <dune/istl/matrix.hh>
16
17
#include <dune/grid/common/grid.hh>

18
#include <dune/amdis/AdaptInfo.hpp>
19
#include <dune/amdis/Assembler.hpp>
20
21
22
23
24
25
26
27
#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>
28
#include <dune/amdis/ProblemStatTraits.hpp>
29
30
#include <dune/amdis/StandardProblemIteration.hpp>

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

37
38
39
#include <dune/amdis/io/FileWriterInterface.hpp>

#include <dune/amdis/utility/DiscreteFunction.hpp>
40
#include <dune/amdis/utility/RangeType.hpp>
41
#include <dune/amdis/utility/TreePath.hpp>
42

43
namespace AMDiS
44
{
45
  template <class GlobalBasis>
46
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
47
      : public ProblemStatBase
48
      , public StandardProblemIteration
49
  {
50
    using Self = ProblemStat;
51

52
  public: // typedefs and static constants
53

54
55
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
56

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

59
    static const int nComponents = 1; // TODO: count leaf nodes in GlobalBasis
60

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

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

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

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

72
73
    using ElementOperator = Operator<GridView, Element>;
    using IntersectionOperator = Operator<GridView, typename GridView::Intersection>;
74

75
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
76

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

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

102
103
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
104
105
    }

106

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

117
118

    /// Adds an operator to \ref A.
119
    /** @{ */
120
121
    template <class RowTreePath, class ColTreePath>
    void addMatrixOperator(ElementOperator const& op, RowTreePath, ColTreePath,
122
123
                           double* factor = nullptr, double* estFactor = nullptr);

124
125
    template <class RowTreePath, class ColTreePath>
    void addMatrixOperator(IntersectionOperator const& op, RowTreePath, ColTreePath,
126
                           double* factor = nullptr, double* estFactor = nullptr);
127

128
    // operator evaluated on the boundary
129
130
    template <class RowTreePath, class ColTreePath>
    void addMatrixOperator(BoundaryType b, IntersectionOperator const& op, RowTreePath, ColTreePath,
131
                           double* factor = nullptr, double* estFactor = nullptr);
132

133
134
135
136
137
138
139
#if 0
    template <class RowTreePath, class ColTreePath>
    void addMatrixOperator(std::map< std::pair<RowTreePath,ColTreePath>, ElementOperator> ops);

    template <class RowTreePath, class ColTreePath>
    void addMatrixOperator(std::map< std::pair<RowTreePath,ColTreePath>, std::pair<ElementOperator, bool> > ops);
#endif
140
    /** @} */
141
142
143


    /// Adds an operator to \ref rhs.
144
    /** @{ */
145
146
    template <class TreePath>
    void addVectorOperator(ElementOperator const& op, TreePath,
147
148
                           double* factor = nullptr, double* estFactor = nullptr);

149
150
    template <class TreePath>
    void addVectorOperator(IntersectionOperator const& op, TreePath,
151
                           double* factor = nullptr, double* estFactor = nullptr);
152

153
    // operator evaluated on the boundary
154
155
    template <class TreePath>
    void addVectorOperator(BoundaryType b, IntersectionOperator const& op, TreePath,
156
                           double* factor = nullptr, double* estFactor = nullptr);
157

158
159
160
161
162
163
164
#if 0
    template <class TreePath>
    void addVectorOperator(std::map<TreePath, ElementOperator> ops);

    template <class TreePath>
    void addVectorOperator(std::map<TreePath, std::pair<ElementOperator, bool> > ops);
#endif
165
    /** @} */
166

167

168
    /// Adds a Dirichlet boundary condition
169
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
170
    void addDirichletBC(Predicate const& predicate,
171
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
172
                        Values const& values);
173

174
175
176
177
178
179
180
181
    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; });
    }

182

183
    /// Implementation of \ref ProblemStatBase::solve
184
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
185
186
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
187
188


189
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
190
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
191
192
193
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
194

195
196

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

Praetorius, Simon's avatar
Praetorius, Simon committed
199

200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
  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; }


241
  public: // get-methods
242

243
    /// Returns a pointer to system-matrix, \ref systemMatrix
244
245
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
246
247


248
    /// Return a pointer to the solution, \ref solution
249
250
    std::shared_ptr<SystemVector>       getSolution()       { return solution; }
    std::shared_ptr<SystemVector const> getSolution() const { return solution; }
251

252
    /// Return a view to a solution component
253
    template <class TreePath>
254
    auto getSolution(TreePath path) const
255
    {
256
      auto tp = makeTreePath(path);
257
      return makeDiscreteFunction(globalBasis, tp, solution);
258
    }
259
260


261
    /// Return a point to the rhs system-vector, \ref rhs
262
263
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
264
265


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

Praetorius, Simon's avatar
Praetorius, Simon committed
269
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
270
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
271
      linearSolver = solver;
272
273
    }

274
275
    /// Return a pointer to the grid, \ref grid
    std::shared_ptr<Grid> getGrid() { return grid; }
276

277
278
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
279
    void setGrid(Grid& grid_)
280
    {
281
282
      grid = std::shared_ptr<Grid>(&grid_, optional_delete(false));
      std::fill(componentGrids.begin(), componentGrids.end(), this->grid.get());
283

284
      createGlobalBasis();
285
286
287
288
289
      createMatricesAndVectors();
      createFileWriter();
    }


290
    /// Return the gridView of the leaf-level
291
    auto leafGridView() { return grid->leafGridView(); }
292
293

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

296
    /// Return the \ref feSpaces
297
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis; }
298
299


300
    /// Implementation of \ref ProblemStatBase::getName
301
302
303
304
    virtual std::string getName() const override
    {
      return name;
    }
305

306
    /// Return a vector of names of the problem components
307
308
309
310
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
311

312
313
314
315
316
    int getNumComponents() const
    {
      return nComponents;
    }

317
  protected: // initialization methods
318

319
    void createGrid()
320
    {
321
322
323
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
      test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
324

325
      grid = MeshCreator<Grid>::create(gridName);
326

327
328
329
330
      msg("Create grid:");
      msg("#elements = "   , grid->size(0));
      msg("#faces/edges = ", grid->size(1));
      msg("#vertices = "   , grid->size(dim));
331
      msg("");
332
    }
333

334
    void createGlobalBasis()
335
    {
336
      globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(leafGridView()));
337
338
339
340
      globalTree = std::make_shared<typename GlobalBasis::LocalView::Tree>(globalBasis->localView().tree());
      matrixOperators.init(*globalTree);
      rhsOperators.init(*globalTree);
      constraints.init(*globalTree);
341
    }
342

343
344
    void createMatricesAndVectors()
    {
345
346
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
      solution = std::make_shared<SystemVector>(*globalBasis, componentNames[0]);
347

348
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
349
      rhs = std::make_shared<SystemVector>(*globalBasis, rhsNames[0]);
350
    }
351

352
    void createSolver()
353
    {
354
355
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
356

357
358
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
359

360
      linearSolver = solverCreator->create(name + "->solver");
361
    }
362

363
    void createFileWriter();
364

Praetorius, Simon's avatar
Praetorius, Simon committed
365

366
  private:
367
368
369
370
371
372
    /// 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;
373

374
375
    /// Grid of this problem.
    std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
376

377
    /// Name of the mesh
378
    std::string gridName = "none";
379

380
    /// Pointer to the meshes for the different problem components
381
    std::vector<Grid*> componentGrids;
382

383
    /// FE spaces of this problem.
384
    std::shared_ptr<GlobalBasis> globalBasis; // eventuell const
385
    std::shared_ptr<typename GlobalBasis::LocalView::Tree> globalTree;
386

387
    /// A FileWriter object
388
    std::list<std::shared_ptr<FileWriterInterface>> filewriter;
389

390
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
391
    std::shared_ptr<LinearSolverType> linearSolver;
392

393
    /// A block-matrix that is filled during assembling
394
    std::shared_ptr<SystemMatrix> systemMatrix;
395

396
    /// A block-vector with the solution components
397
    std::shared_ptr<SystemVector> solution;
398
399

    /// A block-vector (load-vector) corresponding to the right.hand side
400
    /// of the equation, filled during assembling
401
    std::shared_ptr<SystemVector> rhs;
402
403


404
405
  private: // some internal data-structures

406
407
    MatrixOperators<GlobalBasis> matrixOperators;
    VectorOperators<GlobalBasis> rhsOperators;
408

409
    Constraints<GlobalBasis> constraints;
410
  };
411
412


413
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
414
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
415
#endif
416

417
} // end namespace AMDiS
418
419
420

#include "ProblemStat.inc.hpp"