ProblemStat.hpp 14.6 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
38
#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>

#include <amdis/utility/TreePath.hpp>
39

40
namespace AMDiS
41
{
42
  template <class Traits>
43
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
44
      : public ProblemStatBase
45
      , public StandardProblemIteration
46
  {
47
    using Self = ProblemStat;
48

49
  public: // typedefs and static constants
50

51
    using GlobalBasis = typename Traits::GlobalBasis;
52
53
54
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
    using Element  = typename GridView::template Codim<0>::Entity;
55
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
56
57

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

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

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

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

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

78
79
    /// Constructor taking additionally a reference to a mesh that is used
    /// instead of the default created mesh, \ref ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
80
    ProblemStat(std::string name, Grid& grid)
81
      : ProblemStat(std::move(name))
82
    {
83
84
85
86
87
88
89
90
91
92
93
      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_);
94
95
    }

96

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

107
108

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

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


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

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

133

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

140
  public:
141

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


148
149
150
151
152
153
154
155
156
157
158
159
    /// 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;


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

166
167

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

Praetorius, Simon's avatar
Praetorius, Simon committed
170

171
  public: // get-methods
172

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

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

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

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


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


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

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

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

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

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


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

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

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


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

244
  protected: // initialization methods
245

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

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

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

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

    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
269
270
271
272
      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{});
273
    }
274

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

      estimates.resize(leafGridView().indexSet().size(0));
      for (std::size_t i = 0; i < estimates.size(); i++) {
        estimates[i].resize(nComponents);
        for (std::size_t j = 0; j < estimates[i].size(); j++) {
          estimates[i][j] = 0.0; // TODO: Remove when estimate() is implemented
        }
      }
288
    }
289

290
    void createSolver()
291
    {
292
      std::string solverName = "cg";
293
      Parameters::get(name_ + "->solver->name", solverName);
294

295
      auto solverCreator
296
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
297

298
      linearSolver_ = solverCreator->create(name_ + "->solver");
299
    }
300

301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
    void createEstimator()
    {/*
      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));
        }
      }*/
    }


    void createMarker()
    {
      int nMarkersCreated = 0;

      marker.resize(nComponents);

      for (std::size_t i = 0; i < nComponents; i++) {
        marker[i] = Marker<Grid>::
          createMarker(name + "->marker[" + std::to_string(i) + "]", i,
                       estimates, componentGrids[i]);

        if (marker[i]) {
	        nMarkersCreated++;

	        // If there is more than one marker, and all components are defined
	        // on the same mesh, the maximum marking has to be enabled.
 	        if (nMarkersCreated > 1 && nGrids == 1)
 	          marker[i]->setMaximumMarking(true);
        }
      }
    }


348
    void createFileWriter();
349

Praetorius, Simon's avatar
Praetorius, Simon committed
350

351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
  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; }


386
  private:
387
    /// Name of this problem.
388
    std::string name_;
389

390
    /// Grid of this problem.
391
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
392

393
394
395
    /// Number of grids
    int nGrids = 1;

396
    /// Name of the mesh
397
    std::string gridName_ = "none";
398

399
    /// FE spaces of this problem.
400
401
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
402

403
    /// A FileWriter object
404
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
405

406
407
408
409
410
411
    /// Pointer to the adaptation markers
    std::vector<Marker<Grid>* > marker;

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

412
    /// An object of the linearSolver Interface
413
    std::shared_ptr<LinearSolverType> linearSolver_;
414

415
    /// A block-matrix that is filled during assembling
416
    std::shared_ptr<SystemMatrix> systemMatrix_;
417

418
    /// A block-vector with the solution components
419
    std::shared_ptr<SystemVector> solution_;
420

421
422
423
    /// A vector with the local element error estimates
    std::vector<std::vector<double> > estimates;

424
    /// A block-vector (load-vector) corresponding to the right.hand side
425
    /// of the equation, filled during assembling
426
    std::shared_ptr<SystemVector> rhs_;
427
428


429
430
  private: // some internal data-structures

431
432
433
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
434
  };
435

436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
#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};
  }

451

452
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
453
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
454
#endif
455

456
} // end namespace AMDiS
457
458
459

#include "ProblemStat.inc.hpp"