ProblemStat.hpp 13.7 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 Traits>
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
    using GlobalBasis = typename Traits::GlobalBasis;
55
56
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
57

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

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

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

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

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

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

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

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

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

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

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

107

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

118
119

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

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

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

134
135
136
137
138
139
140
#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
141
    /** @} */
142
143
144


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

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

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

159
160
161
162
163
164
165
#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
166
    /** @} */
167

168

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

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

183

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


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

196
197

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

Praetorius, Simon's avatar
Praetorius, Simon committed
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
241
  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; }


242
  public: // get-methods
243

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


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

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


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


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

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

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

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

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


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

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

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


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

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

313
    std::size_t getNumComponents() const
314
315
316
317
    {
      return nComponents;
    }

318
  protected: // initialization methods
319

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

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

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

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

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

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

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

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

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

364
    void createFileWriter();
365

Praetorius, Simon's avatar
Praetorius, Simon committed
366

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

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

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

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

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

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

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

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

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

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


405
406
  private: // some internal data-structures

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

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


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

418
} // end namespace AMDiS
419
420
421

#include "ProblemStat.inc.hpp"