ProblemStat.hpp 13.4 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
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
17
18
#include <dune/grid/common/grid.hh>

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

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

39
40
#include <dune/amdis/utility/RangeType.hpp>

41
namespace AMDiS
42
{
43
  template <class GlobalBasis>
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
53
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
54

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

57
    static const int nComponents = 1; // TODO: count leaf nodes in GlobalBasis
58

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

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

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

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

70
71
    using ElementOperator = Operator<GridView, Element>;
    using IntersectionOperator = Operator<GridView, typename GridView::Intersection>;
72

73
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
74

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

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

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

104

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

115
116

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

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

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

131
132
133
134
135
136
137
#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
138
    /** @} */
139
140
141


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

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

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

156
157
158
159
160
161
162
#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
163
    /** @} */
164

165

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

172

173
    /// Implementation of \ref ProblemStatBase::solve
174
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
175
176
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
177
178


179
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
180
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
181
182
183
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
184

185
186

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

Praetorius, Simon's avatar
Praetorius, Simon committed
189

190
191
192
193
194
195
196
197
198
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
  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; }


231
  public: // get-methods
232

233
    /// Returns a pointer to system-matrix, \ref systemMatrix
234
235
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
236
237


238
    /// Return a pointer to the solution, \ref solution
239
240
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }
241

242
    /// Return the i'th \ref solution component
243
244
    template <class TreePath>
    auto getSolution(TreePath path)
245
    {
246
247
248
249
250
251
252
253
254
255
      using namespace Dune::Functions;

      auto tp = makeTreePath(path);

      using Node = typename GlobalBasis::NodeFactory::template Node<decltype(tp)>;
      using Range = RangeType<Node, double>;

      return makeDiscreteGlobalBasisFunction<Range>(
        subspaceBasis(*globalBasis, tp),
        *solution);
256
    }
257
258


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


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

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

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

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

282
      createGlobalBasis();
283
284
285
286
287
      createMatricesAndVectors();
      createFileWriter();
    }


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

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

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


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

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

310
311
312
313
314
    int getNumComponents() const
    {
      return nComponents;
    }

315
  protected: // initialization methods
316

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

323
      grid = MeshCreator<Grid>::create(gridName);
324

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

332
    void createGlobalBasis()
333
    {
334
335
336
      globalBasis = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(leafGridView()));
      matrixOperators.init(*globalBasis);
      rhsOperators.init(*globalBasis);
337
    }
338

339
340
    void createMatricesAndVectors()
    {
341
342
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
      solution = std::make_shared<SystemVector>(*globalBasis, componentNames[0]);
343

344
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
345
      rhs = std::make_shared<SystemVector>(*globalBasis, rhsNames[0]);
346
    }
347

348
    void createSolver()
349
    {
350
351
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
352

353
354
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
355

356
      linearSolver = solverCreator->create(name + "->solver");
357
    }
358

359
360
    void createFileWriter()
    {
361
      filewriter = std::make_shared<FileWriter<GlobalBasis>>(name + "->output", leafGridView(), componentNames);
362
    }
363

Praetorius, Simon's avatar
Praetorius, Simon committed
364

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

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

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

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

382
    /// FE spaces of this problem.
383
    std::shared_ptr<GlobalBasis> globalBasis; // eventuell const
384

385
    /// A FileWriter object
386
    std::shared_ptr<FileWriter<GlobalBasis>> filewriter;
387

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

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

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

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


402
403
  private: // some internal data-structures

404
405
    MatrixOperators<GlobalBasis> matrixOperators;
    VectorOperators<GlobalBasis> rhsOperators;
406

407
408
    template <class Node>
    struct DirichletData
409
    {
410
      std::list<std::shared_ptr<DirichletBC<WorldVector>>> bc;
411
    };
412

413
    MatrixData<GlobalBasis, DirichletData> matrixConstraints;
414
  };
415
416


417
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
418
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
419
#endif
420

421
} // end namespace AMDiS
422
423
424

#include "ProblemStat.inc.hpp"