ProblemStat.hpp 14.5 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
17
18
#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/DirichletBC.hpp>
19
//#include <amdis/Estimator.hpp>
20
21
22
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
23
#include <amdis/Marker.cpp>
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include <amdis/Mesh.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>

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

#include <amdis/GridFunctions.hpp>
#include <amdis/gridfunctions/DOFVectorView.hpp>

#include <amdis/io/FileWriterInterface.hpp>

38
#include <amdis/utility/TreeData.hpp>
39
#include <amdis/utility/TreePath.hpp>
40

41
namespace AMDiS
42
{
43
  template <class Traits>
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
    using GlobalBasis = typename Traits::GlobalBasis;
53
54
55
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
    using Element  = typename GridView::template Codim<0>::Entity;
56
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
57
58

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

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

64
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
65
    using SystemVector = DOFVector<Traits, double>;
66

67
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
68

69
  public:
70
71
    /**
     * \brief Constructor. Takes the name of the problem that is used to
72
     * access values corresponding to this problem in the parameter file.
73
     **/
74
    explicit ProblemStat(std::string name)
75
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
76
      , name_(std::move(name))
77
    {}
78

79
80
    /// Constructor taking additionally a reference to a grid that is used
    /// instead of the default created grid, \ref ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
81
    ProblemStat(std::string name, Grid& grid)
82
      : ProblemStat(std::move(name))
83
    {
84
85
86
87
88
89
90
91
92
93
94
      grid_ = Dune::stackobject_to_shared_ptr(grid);
      Parameters::get(name_ + "->mesh", gridName_);
    }

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
    ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
      : ProblemStat(std::move(name), grid)
    {
      globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
      initGlobalBasis(*globalBasis_);
95
96
    }

97

98
    /**
99
     * \brief Initialisation of the problem.
100
     *
101
102
103
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
104
    void initialize(Flag initFlag,
105
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
106
                    Flag adoptFlag = INIT_NOTHING);
107

108
109

    /// Adds an operator to \ref A.
110
    /** @{ */
111
112
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath, ColTreePath,
113
                           double* factor = nullptr, double* estFactor = nullptr);
114

115
    // operator evaluated on the boundary
116
117
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath, ColTreePath,
118
                           double* factor = nullptr, double* estFactor = nullptr);
119
    /** @} */
120
121
122


    /// Adds an operator to \ref rhs.
123
    /** @{ */
124
125
    template <class Operator, class TreePath>
    void addVectorOperator(Operator const& op, TreePath,
126
                           double* factor = nullptr, double* estFactor = nullptr);
127

128
    // operator evaluated on the boundary
129
130
    template <class Operator, class TreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath,
131
                           double* factor = nullptr, double* estFactor = nullptr);
132
    /** @} */
133

134

135
    /// Adds a Dirichlet boundary condition
136
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
137
    void addDirichletBC(Predicate const& predicate,
138
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
139
                        Values const& values);
140

141
  public:
142

143
    /// Implementation of \ref ProblemStatBase::solve
144
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
145
146
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
147
148


149
150
151
152
153
154
155
156
157
158
159
160
    /// Implementation of \ref ProblemStatBase::estimate.
    virtual void estimate(AdaptInfo& adaptInfo) override;


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


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


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

167
168

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

Praetorius, Simon's avatar
Praetorius, Simon committed
171

172
  public: // get-methods
173

174
175
176
    /// Returns a pointer to system-matrix, \ref systemMatrix_
    std::shared_ptr<SystemMatrix> getSystemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
177

178
    /// Returns a pointer to the solution vector, \ref solution_
179
    std::shared_ptr<SystemVector> getSolutionVector() const
180
    {
181
      return solution_;
182
183
184
185
186
187
188
    }

    /// Return a mutable view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {})
    {
      auto&& tp = makeTreePath(path);
189
      return makeDOFVectorView(*solution_, tp);
190
    }
191

192
193
194
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
195
    {
196
      auto&& tp = makeTreePath(path);
197
      return makeDOFVectorView(*solution_, tp);
198
    }
199
200


201
    /// Return a point to the rhs system-vector, \ref rhs
202
203
    std::shared_ptr<SystemVector> getRhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
204
205


206
    /// Return a pointer to the linear solver, \ref linearSolver
207
    std::shared_ptr<LinearSolverType> getSolver() { return linearSolver_; }
208

Praetorius, Simon's avatar
Praetorius, Simon committed
209
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
210
    {
211
      linearSolver_ = solver;
212
213
    }

214
    /// Return a pointer to the grid, \ref grid
215
    std::shared_ptr<Grid> getGrid() { return grid_; }
216

217
    /// Set the grid. Stores pointer to passed reference and initializes feSpaces
218
    /// matrices and vectors, as well as the file-writer.
219
    void setGrid(Grid& grid)
220
    {
221
222
      grid_ = Dune::stackobject_to_shared_ptr(grid);

223
      createGlobalBasis();
224
225
226
227
228
      createMatricesAndVectors();
      createFileWriter();
    }


229
    /// Return the gridView of the leaf-level
230
231
232
233
    auto leafGridView() { return grid_->leafGridView(); }

    /// Return the gridView of levle `level`
    auto levelGridView(int level) { return grid_->levelGridView(level); }
234

235
    /// Return the \ref feSpaces
236
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
237
238


239
    /// Implementation of \ref ProblemStatBase::getName
240
241
    virtual std::string getName() const override
    {
242
      return name_;
243
    }
244

245
  protected: // initialization methods
246

247
    void createGrid()
248
    {
249
250
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
251
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
252

253
      grid_ = MeshCreator<Grid>::create(gridName_);
254

255
      msg("Create grid:");
256
257
258
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
259
      msg("");
260
    }
261

262
    void createGlobalBasis()
263
    {
264
265
      globalBasis_ = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid_->leafGridView()));
      initGlobalBasis(*globalBasis_);
266
267
268
269
    }

    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
270
271
272
273
      localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
      matrixOperators_.init(localView_->tree(), tag::store{});
      rhsOperators_.init(localView_->tree(), tag::store{});
      constraints_.init(localView_->tree(), tag::store{});
274
    }
275

276
277
    void createMatricesAndVectors()
    {
278
279
280
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
281

282
283
284
285
286
287
288
289
290
      systemMatrix = std::make_shared<SystemMatrix>(*globalBasis, *globalBasis, "mat");
      solution = std::make_shared<SystemVector>(*globalBasis, "solution");

      auto localView = globalBasis->localView();
      forEachNode(localView.tree(), [&,this](auto const& node, auto treePath)
      {
        std::string i = to_string(treePath);
        estimates[i].resize(globalBasis->gridView().indexSet().size(0));
        for (std::size_t j = 0; j < estimates[i].size(); j++)
291
          estimates[i][j] = 0.0; // TODO: Remove when estimate() is implemented
292
293
294
      });

      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
295
    }
296

297
    void createSolver()
298
    {
299
      std::string solverName = "cg";
300
      Parameters::get(name_ + "->solver->name", solverName);
301

302
      auto solverCreator
303
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
304

305
      linearSolver_ = solverCreator->create(name_ + "->solver");
306
    }
307

308
    void createEstimator()
309
    {/* COPIED FROM AMDiS
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
      for (std::size_t i = 0, i < nComponents, i++) {
        std::string estName = "";
        Parameters::get(name + "->estimator->name[" + std::to_string(i) + "]", estName);
        auto estimatorCreator
          = named(CreatorMap<Estimator>::getCreator(estName, name + "->estimator->name[" + std::to_string(i) + "]"));
        if (estimatorCreator) {
	        estimatorCreator->setName(estName);
	        estimatorCreator->setRow(i);
	        estimatorCreator->setSolution(solution->getDOFVector(i));
	        estimator[i] = estimatorCreator->create();
        }

        if (estimator[i]) {
          for (std::size_t j = 0; j < nComponents; j++)
	          estimator[i]->addSystem((*systemMatrix)[i][j],
              solution->getDOFVector(j),
				      rhs->getDOFVector(i));
        }
      }*/
    }


332
    void createMarker();
333

334
    void createFileWriter();
335

Praetorius, Simon's avatar
Praetorius, Simon committed
336

337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
  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
    {
      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; }


372
  private:
373
    /// Name of this problem.
374
    std::string name_;
375

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

379
380
381
    /// Number of grids
    int nGrids = 1;

382
    /// Name of the grid
383
    std::string gridName_ = "none";
384

385
    /// FE spaces of this problem.
386
387
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
388

389
    /// A FileWriter object
390
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
391

392
    /// Pointer to the adaptation markers
393
    std::list<Marker<Traits>* > marker;
394
395
396
397

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

398
    /// An object of the linearSolver Interface
399
    std::shared_ptr<LinearSolverType> linearSolver_;
400

401
    /// A block-matrix that is filled during assembling
402
    std::shared_ptr<SystemMatrix> systemMatrix_;
403

404
    /// A block-vector with the solution components
405
    std::shared_ptr<SystemVector> solution_;
406

407
    /// A vector with the local element error estimates
408
    /// reverse indexed by [component index][element index]
409
    std::map<std::string, std::vector<double> > estimates;
410

411
    /// A block-vector (load-vector) corresponding to the right.hand side
412
    /// of the equation, filled during assembling
413
    std::shared_ptr<SystemVector> rhs_;
414
415


416
417
  private: // some internal data-structures

418
419
420
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
421
422

//    VectorData<GlobalBasis, TODO: NodeData> estimates_;
423
  };
424

425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
  // Deduction rule
  template <class Grid, class GlobalBasis>
  ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
    -> ProblemStat<DefaultProblemTraits<GlobalBasis>>;
#endif

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

440

441
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
442
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
443
#endif
444

445
} // end namespace AMDiS
446
447
448

#include "ProblemStat.inc.hpp"