ProblemStat.hpp 13.8 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
#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
26
#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>
27
#include <dune/amdis/ProblemStatTraits.hpp>
28
29
#include <dune/amdis/StandardProblemIteration.hpp>

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

36
namespace AMDiS
37
{
38
  struct BoundaryType { int b; };
39

40
  template <class GlobalBasis>
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
50
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
51

52
    using Element  = typename GridView::template Codim<0>::Entity;
53

54

55
    /// Dimension of the mesh
56
    static constexpr int dim = Grid::dimension;
57

58
    /// Dimension of the world
59
    static constexpr int dow = Grid::dimensionworld;
60

61
    /// Number of problem components
62
    static constexpr int nComponents = Traits::nComponents; /* number of leaf nodes in the basis-tree */
63
64


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

67
68
    using SystemVectorType = SystemVector<GlobalBasis>;
    using SystemMatrixType = SystemMatrix<GlobalBasis>;
69

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

73
74
    using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
                                                   typename SystemVectorType::MultiVector>;
75

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

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

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

105

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

116
117

    /// Adds an operator to \ref A.
118
    /** @{ */
119
120
121
122
    void addMatrixOperator(ElementOperator const& op, int i, int j,
                           double* factor = nullptr, double* estFactor = nullptr);

    void addMatrixOperator(IntersectionOperator const& op, int i, int j,
123
                           double* factor = nullptr, double* estFactor = nullptr);
124

125
    // operator evaluated on the boundary
126
    void addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j,
127
                           double* factor = nullptr, double* estFactor = nullptr);
128

Praetorius, Simon's avatar
Praetorius, Simon committed
129
130
    void addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops);
    void addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator, bool> > ops);
131
    /** @} */
132
133
134


    /// Adds an operator to \ref rhs.
135
    /** @{ */
136
137
138
139
    void addVectorOperator(ElementOperator const& op, int i,
                           double* factor = nullptr, double* estFactor = nullptr);

    void addVectorOperator(IntersectionOperator const& op, int i,
140
                           double* factor = nullptr, double* estFactor = nullptr);
141

142
    // operator evaluated on the boundary
143
    void addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i,
144
                           double* factor = nullptr, double* estFactor = nullptr);
145

Praetorius, Simon's avatar
Praetorius, Simon committed
146
147
    void addVectorOperator(std::map< int, ElementOperator> ops);
    void addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops);
148
    /** @} */
149

150

151
152
    /// Adds a Dirichlet boundary condition
    template <class Predicate, class Values>
153
    void addDirichletBC(Predicate const& predicate,
Praetorius, Simon's avatar
Praetorius, Simon committed
154
155
                        int row, int col,
                        Values const& values);
156

157

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


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

170
171

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

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


216
  public: // get-methods
217
218

    /// Returns nr of components \ref nComponents
219
    static constexpr int getNumComponents()
220
221
222
    {
      return nComponents;
    }
223
224


225
    /// Returns a pointer to system-matrix, \ref systemMatrix
226
227
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }
228
229


230
    /// Return a pointer to the solution, \ref solution
231
232
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }
233

234
    /// Return the i'th \ref solution component
235
236
    template <std::size_t I = 0>
    decltype(auto) getSolution(const index_t<I> _i = {})
237
    {
238
      return (*solution)[_i];
239
    }
240
241


242
    /// Return a point to the rhs system-vector, \ref rhs
243
244
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }
245
246


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

Praetorius, Simon's avatar
Praetorius, Simon committed
250
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
251
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
252
      linearSolver = solver;
253
254
    }

255
256
    /// Return a pointer to the grid, \ref grid
    std::shared_ptr<Grid> getGrid() { return grid; }
257

258
259
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
260
    void setGrid(Grid& grid_)
261
    {
262
263
      grid = std::shared_ptr<Grid>(&grid_, optional_delete(false));
      std::fill(componentGrids.begin(), componentGrids.end(), this->grid.get());
264
265
266
267
268
269
270

      createFeSpaces();
      createMatricesAndVectors();
      createFileWriter();
    }


271
    /// Return the gridView of the leaf-level
272
    auto leafGridView() { return grid->leafGridView(); }
273
274

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

277
    /// Return the \ref feSpaces
278
    auto getGlobalBasis() { return globalBasis; }
279
280


281
    /// Return the I'th \ref feSpaces component
282
    template <std::size_t I = 0>
283
    FeSpace<I> const& getComponentBasis(const index_t<I> = {}) const
284
    {
285
      return std::get<I>(*globalBasis);
286
287
288
    }

    template <std::size_t I = 0>
289
    FeSpace<I>& getComponentBasis(const index_t<I> = {})
290
    {
291
      return std::get<I>(*globalBasis);
292
    }
293

294
    /// Implementation of \ref ProblemStatBase::getName
295
296
297
298
    virtual std::string getName() const override
    {
      return name;
    }
299

300
    /// Return a vector of names of the problem components
301
302
303
304
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }
305

306
  protected: // initialization methods
307

308
    void createGrid()
309
    {
310
311
312
      gridName = "";
      Parameters::get(name + "->mesh", gridName);
      test_exit(!gridName.empty(), "No mesh name specified for '", name, "->mesh'!");
313

314
      grid = MeshCreator<Grid>::create(gridName);
315

316
317
318
319
      msg("Create grid:");
      msg("#elements = "   , grid->size(0));
      msg("#faces/edges = ", grid->size(1));
      msg("#vertices = "   , grid->size(dim));
320
      msg("");
321
    }
322

323
    void createGlobalBasis()
324
    {
325
      globalBasis = std::make_shared<GlobalBasis>(new GlobalBasis(leafGridView()));
326
    }
327

328
329
    void createMatricesAndVectors()
    {
330
331
      systemMatrix = std::make_shared<SystemMatrixType>(*globalBasis);
      solution = std::make_shared<SystemVectorType>(*globalBasis, componentNames);
332

333
      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
334
      rhs = std::make_shared<SystemVectorType>(*globalBasis, rhsNames);
335
    }
336

337
    void createSolver()
338
    {
339
340
      std::string solverName = "cg";
      Parameters::get(name + "->solver->name", solverName);
341

342
343
      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name + "->solver->name"));
344

345
      linearSolver = solverCreator->create(name + "->solver");
346
    }
347

348
349
    void createFileWriter()
    {
350
      filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", leafGridView(), componentNames);
351
    }
352

Praetorius, Simon's avatar
Praetorius, Simon committed
353

354
  private:
355
356
357
358
359
360
    /// 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;
361

362
363
    /// Grid of this problem.
    std::shared_ptr<Grid> grid; // TODO: generalize to multi-mesh problems
364

365
    /// Name of the mesh
366
    std::string gridName = "none";
367

368
    /// Pointer to the meshes for the different problem components
369
    std::vector<Grid*> componentGrids;
370

371
    /// FE spaces of this problem.
372
    std::shared_ptr<GlobalBasis> globalBasis; // eventuell const
373

374
    /// A FileWriter object
Praetorius, Simon's avatar
Praetorius, Simon committed
375
    std::shared_ptr<FileWriter<Traits>> filewriter;
376

377
    /// An object of the linearSolver Interface
Praetorius, Simon's avatar
Praetorius, Simon committed
378
    std::shared_ptr<LinearSolverType> linearSolver;
379

380
    /// A block-matrix that is filled during assembling
Praetorius, Simon's avatar
Praetorius, Simon committed
381
    std::shared_ptr<SystemMatrixType> systemMatrix;
382

383
    /// A block-vector with the solution components
Praetorius, Simon's avatar
Praetorius, Simon committed
384
    std::shared_ptr<SystemVectorType> solution;
385
386

    /// A block-vector (load-vector) corresponding to the right.hand side
387
    /// of the equation, filled during assembling
Praetorius, Simon's avatar
Praetorius, Simon committed
388
    std::shared_ptr<SystemVectorType> rhs;
389
390


391
392
  private: // some internal data-structures

393
    template <class T>
394
    using MatrixEntries = Dune::FieldMatrix<T, nComponents, nComponents>;
395

396
    template <class T>
397
    using VectorEntries = Dune::FieldVector<T, nComponents>;
398

399
400
401
402
403

    /// Helper class to store an operator with optionally factors and
    /// a boundary indicator
    template <class OperatorType>
    struct Scaled
404
    {
405
406
407
408
      std::shared_ptr<OperatorType> op;
      double* factor = nullptr;
      double* estFactor = nullptr;
      BoundaryType b = {0};
409
    };
410

411
412
    // lists of operators on elements, boundary and intersection
    struct Operators
413
    {
414
415
416
      std::list<Scaled<ElementOperator>> element;
      std::list<Scaled<IntersectionOperator>> boundary;
      std::list<Scaled<IntersectionOperator>> intersection;
417

418
419
420
421
      std::list<std::shared_ptr<DirichletBC<WorldVector>>> dirichlet; // boundary conditions

      bool assembled = false; // if false, do reassemble
      bool changing = false; // if true, or vectorAssembled false, do reassemble
422
    };
423
424
425

    MatrixEntries<Operators> matrixOperators;
    VectorEntries<Operators> rhsOperators;
426
  };
427
428


429
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
430
431
  extern template class ProblemStat<TestTraits<2,1>>; // 1 component with lagrange degree 1
  extern template class ProblemStat<TestTraits<2,1,2>>; // 2 components with different polynomial degree
432
#endif
433

434
} // end namespace AMDiS
435
436
437

#include "ProblemStat.inc.hpp"