ProblemStat.hpp 17.2 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
14
#include <dune/grid/common/grid.hh>

15
16
#include <amdis/AdaptInfo.hpp>
#include <amdis/CreatorInterface.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
17
#include <amdis/CreatorMap.hpp>
18
#include <amdis/DirichletBC.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
19
//#include <amdis/Estimator.hpp>
20
21
22
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
23
#include <amdis/LinearSolvers.hpp>
24
#include <amdis/LocalAssemblerList.hpp>
25
#include <amdis/Marker.hpp>
26
#include <amdis/Mesh.hpp>
27
#include <amdis/PeriodicBC.hpp>
28
29
30
31
32
33
34
35
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>

#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/Utility.hpp>

#include <amdis/GridFunctions.hpp>
36
#include <amdis/gridfunctions/DiscreteFunction.hpp>
37
38
39
40
#include <amdis/gridfunctions/DOFVectorView.hpp>

#include <amdis/io/FileWriterInterface.hpp>

41
#include <amdis/utility/TreeData.hpp>
42
#include <amdis/utility/TreePath.hpp>
43

44
namespace AMDiS
45
{
46
47
48
49
  // forward declaration
  template <class Traits>
  class ProblemInstat;

50
  template <class Traits>
51
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
52
      : public ProblemStatBase
53
      , public StandardProblemIteration
54
  {
55
    using Self = ProblemStat;
56

57
58
    friend class ProblemInstat<Traits>;

59
  public: // typedefs and static constants
60

61
    using GlobalBasis = typename Traits::GlobalBasis;
Praetorius, Simon's avatar
Praetorius, Simon committed
62
63
64
    using GridView    = typename GlobalBasis::GridView;
    using Grid        = typename GridView::Grid;
    using Element     = typename GridView::template Codim<0>::Entity;
65
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
66
    using WorldMatrix = FieldMatrix<typename WorldVector::field_type, WorldVector::dimension, WorldVector::dimension>;
67

68
    /// Dimension of the grid
69
    static constexpr int dim = Grid::dimension;
70

71
    /// Dimension of the world
72
    static constexpr int dow = Grid::dimensionworld;
73

74
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
75
    using SystemVector = DOFVector<GlobalBasis, double>;
76

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

79
  public:
80
81
    /**
     * \brief Constructor. Takes the name of the problem that is used to
82
     * access values corresponding to this problem in the parameter file.
83
     **/
84
    explicit ProblemStat(std::string const& name)
85
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
86
      , name_(name)
87
    {}
88

89
90
    /// Constructor taking additionally a reference to a grid that is used
    /// instead of the default created grid, \ref ProblemStat
91
92
    ProblemStat(std::string const& name, Grid& grid)
      : ProblemStat(name)
93
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
94
      adoptGrid(Dune::stackobject_to_shared_ptr(grid));
95
96
97
98
    }

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
99
100
    ProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
      : ProblemStat(name, grid)
101
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
102
      adoptGlobalBasis(Dune::stackobject_to_shared_ptr(globalBasis));
103
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
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
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
112
    void initialize(Flag initFlag, Self* adoptProblem = nullptr, Flag adoptFlag = INIT_NOTHING);
113

114

Praetorius, Simon's avatar
Praetorius, Simon committed
115
    /// Add an operator to \ref A.
116
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
    /// Operator evaluated on the whole element
    /**
     * Adds an operator to the list of element operators to be assembled in
     * quadrature points inside the element.
     *
     * \param op   A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator
     * \param row  TreePath identifying the sub-basis in the global basis tree
     *             corresponding to the row basis. \see treepath()
     * \param col  TreePath identifying the sub-basis in the global basis tree
     *             corresponding to the column basis. \see treepath()
     *
     * Example:
     * ```
     * auto op = makeOperator(tag::test_trial{}, 1.0/tau);
     * prob.addMatrixOperator(op, _0, _0);
     * ```
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
134
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
135
136
137
138
    void addMatrixOperator(Operator const& op, RowTreePath row = {}, ColTreePath col = {})
    {
      systemMatrix_->addOperator(tag::element_operator<Element>{}, op, row, col);
    }
139

Praetorius, Simon's avatar
Praetorius, Simon committed
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
    /// Operator evaluated on the boundary of the domain with boundary index `b`
    /**
     * Adds an operator to the list of boundary operators to be assembled in
     * quadrature points on the boundary intersections.
     *
     * \param b    Boundary indentifier where to assemble this operator. Can be
     *             constructed from an integer. \see BoundaryType
     * \param op   A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator
     * \param row  TreePath identifying the sub-basis in the global basis tree
     *             corresponding to the row basis. \see treepath()
     * \param col  TreePath identifying the sub-basis in the global basis tree
     *             corresponding to the column basis. \see treepath()
     *
     * Example:
     * ```
     * auto op = makeOperator(tag::test_trial{}, alpha);
     * prob.addMatrixOperator(BoundaryType{1}, op, _0, _0);
     * ```
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
159
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
160
161
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath row = {}, ColTreePath col = {})
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
162
      using I = typename GridView::Intersection;
163
      systemMatrix_->addOperator(tag::boundary_operator<I>{boundaryManager_,b}, op, row, col);
164
    }
165
    /** @} */
166
167


Praetorius, Simon's avatar
Praetorius, Simon committed
168
    /// Add an operator to \ref rhs.
169
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
    /// Operator evaluated on the whole element
    /**
     * Adds an operator to the list of element operators to be assembled in
     * quadrature points inside the element.
     *
     * \param op    A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator
     * \param path  TreePath identifying the sub-basis in the global basis tree
     *              corresponding to the row basis. \see treepath()
     *
     * Example:
     * ```
     * auto op = makeOperator(tag::test{}, probInstat.getOldSolution(0) / tau);
     * prob.addVectorOperator(op, _0);
     * ```
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
185
    template <class Operator, class TreePath = RootTreePath>
186
187
188
189
    void addVectorOperator(Operator const& op, TreePath path = {})
    {
      rhs_->addOperator(tag::element_operator<Element>{}, op, path);
    }
190

Praetorius, Simon's avatar
Praetorius, Simon committed
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
    /// Operator evaluated on the boundary of the domain with boundary index `b`
    /**
     * Adds an operator to the list of boundary operators to be assembled in
     * quadrature points on the boundary intersections.
     *
     * \param b     Boundary indentifier where to assemble this operator. Can be
     *              constructed from an integer. \see BoundaryType
     * \param op    A (pre-) local operator, \see LocalOperator, \see GridFunctionOperator
     * \param path  TreePath identifying the sub-basis in the global basis tree
     *              corresponding to the row basis. \see treepath()
     *
     * Example:
     * ```
     * auto op = makeOperator(tag::test{}, [g](auto const& x) { return g(x); });
     * prob.addVectorOperator(BoundaryType{1}, op, _0);
     * ```
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
208
    template <class Operator, class TreePath = RootTreePath>
209
210
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath path = {})
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
211
      using I = typename GridView::Intersection;
212
      rhs_->addOperator(tag::boundary_operator<I>{boundaryManager_,b}, op, path);
213
    }
214
    /** @} */
215

216

Praetorius, Simon's avatar
Praetorius, Simon committed
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
    /// Add boundary conditions to the system
    /** @{ */
    /// Dirichlet boundary condition
    /**
     * Enforce Dirichlet boundary values for the solution vector on boundary
     * regions identified by the predicate.
     *
     * \param predicate  Functor `bool(WorldVector)` returning true for all
     *                   DOFs on the boundary that should be assigned a value.
     * \param row        TreePath identifying the sub-basis in the global basis tree
     *                   corresponding to the row basis. \see treepath()
     * \param col        TreePath identifying the sub-basis in the global basis tree
     *                   corresponding to the column basis. \see treepath()
     * \param values     Functor `Range(WorldVector)` or any \ref GridFunction
     *                   that is evaluated in the DOFs identified by the predicate.
     *
     * Example:
     * ```
     * prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8; }, 0, 0,
     *                     [](auto const& x) { return 0.0; });
     * ```
     **/
239
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
240
    void addDirichletBC(Predicate const& predicate,
241
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
242
                        Values const& values);
243
244
245
246
247
248
249
250
251

    template <class RowTreePath, class ColTreePath, class Values>
    void addDirichletBC(BoundaryType id,
                        RowTreePath row, ColTreePath col,
                        Values const& values);

    /// Add a periodic boundary conditions to the system, by specifying a face transformation
    /// y = A*x + b of coordinates. We assume, that A is orthonormal.
    void addPeriodicBC(BoundaryType id, WorldMatrix const& A, WorldVector const& b);
Praetorius, Simon's avatar
Praetorius, Simon committed
252
253
    /** @} */

254

255
  public:
256

257
    /// Implementation of \ref ProblemStatBase::solve
258
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
259
260
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
261

262
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
263
264
265
266
    virtual void buildAfterAdapt(AdaptInfo& adaptInfo,
                                 Flag flag,
                                 bool asmMatrix = true,
                                 bool asmVector = true) override;
267

Praetorius, Simon's avatar
Praetorius, Simon committed
268
269
270
271
272
273
274
    /// \brief Assemble the linear system by calling \ref buildAfterAdapt with
    /// `asmMatrix` and `asmVector` set to true.
    void assemble(AdaptInfo& adaptInfo)
    {
      buildAfterAdapt(adaptInfo, Flag{0}, true, true);
    }

275
    /// Writes output files.
276
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
277

Praetorius, Simon's avatar
Praetorius, Simon committed
278

279
  public: // get-methods
280

Praetorius, Simon's avatar
Praetorius, Simon committed
281
282
283
284
285
286
287
288
289
    /// Implementation of \ref ProblemStatBase::name
    virtual std::string const& name() const override { return name_; }


    /// Return a reference to the grid, \ref grid
    Grid&       grid()       { return *grid_; }
    Grid const& grid() const { return *grid_; }

    /// Return the gridView of the leaf-level
290
    GridView const& gridView() const { return globalBasis_->gridView(); }
Praetorius, Simon's avatar
Praetorius, Simon committed
291

292
293
294
295
296
297
    /// Return the boundary manager to identify boundary segments
    BoundaryManager<Grid>&       boundaryManager()       { return *boundaryManager_; }
    BoundaryManager<Grid> const& boundaryManager() const { return *boundaryManager_; }

    auto const& boundaryManagerPtr() { return boundaryManager_; }

Praetorius, Simon's avatar
Praetorius, Simon committed
298
299
300
301
302
303
304
305
    /// Return the \ref globalBasis_
    GlobalBasis&       globalBasis()       { return *globalBasis_; }
    GlobalBasis const& globalBasis() const { return *globalBasis_; }

    /// Return a reference to the linear solver, \ref linearSolver
    LinearSolverType&       solver()       { return *linearSolver_; }
    LinearSolverType const& solver() const { return *linearSolver_; }

306
    /// Returns a reference to system-matrix, \ref systemMatrix_
Praetorius, Simon's avatar
Praetorius, Simon committed
307
308
    SystemMatrix&       systemMatrix()       { return *systemMatrix_; }
    SystemMatrix const& systemMatrix() const { return *systemMatrix_; }
309
310

    /// Returns a reference to the solution vector, \ref solution_
Praetorius, Simon's avatar
Praetorius, Simon committed
311
312
    SystemVector&       solutionVector()       { return *solution_; }
    SystemVector const& solutionVector() const { return *solution_; }
313
314

    /// Return a reference to the rhs system-vector, \ref rhs
Praetorius, Simon's avatar
Praetorius, Simon committed
315
316
    SystemVector&       rhsVector()       { return *rhs_; }
    SystemVector const& rhsVector() const { return *rhs_; }
317

318
319
320

    /// Return a mutable view to a solution component
    template <class TreePath = RootTreePath>
Praetorius, Simon's avatar
Praetorius, Simon committed
321
    auto solution(TreePath path = {})
322
323
    {
      auto&& tp = makeTreePath(path);
324
      return makeDOFVectorView(*solution_, tp);
325
    }
326

327
328
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
Praetorius, Simon's avatar
Praetorius, Simon committed
329
    auto solution(TreePath path = {}) const
330
    {
331
      auto&& tp = makeTreePath(path);
332
      return makeDiscreteFunction(*solution_, tp);
333
    }
334
335


Praetorius, Simon's avatar
Praetorius, Simon committed
336
  public: // set-methods
337

338
    /// Set a new linear solver for the problem
Praetorius, Simon's avatar
Praetorius, Simon committed
339
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
340
    {
341
      linearSolver_ = solver;
342
343
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
344
345
346
347
348
    void setSolver(LinearSolverType& solver)
    {
      setSolver(Dune::stackobject_to_shared_ptr(solver));
    }

349

350
    /// Set the grid. Stores pointer and initializes feSpaces
Praetorius, Simon's avatar
Praetorius, Simon committed
351
    /// matrices and vectors, as well as markers and file-writers.
352
    void setGrid(std::shared_ptr<Grid> const& grid)
353
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
354
      adoptGrid(grid);
355
      createGlobalBasis();
356
      createMatricesAndVectors();
Praetorius, Simon's avatar
Praetorius, Simon committed
357
      createMarker();
358
359
360
      createFileWriter();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
361
362
363
364
365
    void setGrid(Grid& grid)
    {
      setGrid(Dune::stackobject_to_shared_ptr(grid));
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
366
367
368
369
370
371
372
373
374
375
376
377
    void addMarker(std::shared_ptr<Marker<Grid>> const& marker)
    {
      marker_.push_back(marker);
      if (marker_.size() > 1)
        marker_.back()->setMaximumMarking(true);
    }

    void addMarker(Marker<Grid>& marker)
    {
      addMarker(Dune::stackobject_to_shared_ptr(marker));
    }

378

379
  protected: // initialization methods
380

Praetorius, Simon's avatar
Praetorius, Simon committed
381
382
383
384
385
386
    void createGlobalBasis();
    void createGrid();
    void createMatricesAndVectors();
    void createSolver();
    void createMarker();
    void createFileWriter();
387

Praetorius, Simon's avatar
Praetorius, Simon committed
388
    void adoptGlobalBasis(std::shared_ptr<GlobalBasis> const& globalBasis)
389
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
390
      globalBasis_ = globalBasis;
391
392
393
      initGlobalBasis(*globalBasis_);
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
394
    void adoptGrid(std::shared_ptr<Grid> const& grid)
395
396
397
398
399
400
    {
      adoptGrid(grid, std::make_shared<BoundaryManager<Grid>>(grid));
    }

    void adoptGrid(std::shared_ptr<Grid> const& grid,
                   std::shared_ptr<BoundaryManager<Grid>> const& boundaryManager)
401
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
402
      grid_ = grid;
403
      boundaryManager_ = boundaryManager;
404
      Parameters::get(name_ + "->mesh", gridName_);
405
    }
406

Praetorius, Simon's avatar
Praetorius, Simon committed
407
  private:
408

Praetorius, Simon's avatar
Praetorius, Simon committed
409
410
    void createGlobalBasisImpl(std::true_type);
    void createGlobalBasisImpl(std::false_type);
411

Praetorius, Simon's avatar
Praetorius, Simon committed
412
    void initGlobalBasis(GlobalBasis const& globalBasis);
413

Praetorius, Simon's avatar
Praetorius, Simon committed
414

415
416
  public: // implementation of iteration interface methods

Praetorius, Simon's avatar
Praetorius, Simon committed
417
    /// Implementation of \ref StandardProblemIteration::oneIteration.
418
419
420
421
422
423
    virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
    {
      return StandardProblemIteration::oneIteration(adaptInfo, toDo);
    }

    /// Implementation of \ref ProblemStatBase::estimate.
Praetorius, Simon's avatar
Praetorius, Simon committed
424
    virtual void estimate(AdaptInfo& adaptInfo) override { /* do nothing. */ }
425
426

    /// Implementation of \ref ProblemStatBase::refineMesh.
427
    virtual Flag adaptGrid(AdaptInfo& adaptInfo) override;
428
429

    /// Implementation of \ref ProblemStatBase::markElements.
430
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
431

432
433
434
435
436
    /// Uniform global grid coarsening by up to n level
    virtual Flag globalCoarsen(int n) override;

    /// Uniform global refinement by n level
    virtual Flag globalRefine(int n) override;
437

438
  private:
439
    /// Name of this problem.
440
    std::string name_;
441

442
    /// Grid of this problem.
Praetorius, Simon's avatar
Praetorius, Simon committed
443
    std::shared_ptr<Grid> grid_;
444

445
    /// Name of the grid
446
    std::string gridName_ = "mesh";
447

448
449
450
    /// Management of boundary conditions
    std::shared_ptr<BoundaryManager<Grid>> boundaryManager_;

451
    /// FE spaces of this problem.
452
    std::shared_ptr<GlobalBasis> globalBasis_;
453

454
    /// A FileWriter object
455
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
456

457
    /// Pointer to the adaptation markers
Praetorius, Simon's avatar
Praetorius, Simon committed
458
    std::list<std::shared_ptr<Marker<Grid>>> marker_;
Praetorius, Simon's avatar
Praetorius, Simon committed
459
460
461

    /// Pointer to the estimators for this problem
//    std::vector<Estimator*> estimator;
462

463
    /// An object of the linearSolver Interface
464
    std::shared_ptr<LinearSolverType> linearSolver_;
465

Praetorius, Simon's avatar
Praetorius, Simon committed
466
    /// Matrix that is filled during assembling
467
    std::shared_ptr<SystemMatrix> systemMatrix_;
468

Praetorius, Simon's avatar
Praetorius, Simon committed
469
    /// Vector with the solution components
470
    std::shared_ptr<SystemVector> solution_;
471

Praetorius, Simon's avatar
Praetorius, Simon committed
472
    /// Vector (load-vector) corresponding to the right-hand side
473
    /// of the equation, filled during assembling
474
    std::shared_ptr<SystemVector> rhs_;
475

Praetorius, Simon's avatar
Praetorius, Simon committed
476
477
478
479
    /// A vector with the local element error estimates
    /// for each node in the basis tree, indexed by [to_string(treePath)][element index]
    std::map<std::string, std::vector<double>> estimates_;

480

481
482
  private: // some internal data-structures

483
484
    DirichletBCs<GlobalBasis, GlobalBasis> dirichletBCs_;
    PeriodicBCs<GlobalBasis, GlobalBasis> periodicBCs_;
485
  };
486

Praetorius, Simon's avatar
Praetorius, Simon committed
487

488
489
490
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
  // Deduction rule
  template <class Grid, class GlobalBasis>
491
  ProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
492
493
494
495
496
497
    -> ProblemStat<DefaultProblemTraits<GlobalBasis>>;
#endif

  // Generator for ProblemStat with given Grid and GlobalBasis
  template <class Grid, class GlobalBasis>
  ProblemStat<DefaultProblemTraits<GlobalBasis>>
498
  makeProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
499
  {
500
    return {name, grid, globalBasis};
501
502
  }

503
} // end namespace AMDiS
504
505
506

#include "ProblemStat.inc.hpp"