ProblemStat.hpp 13.1 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 19 20 21
#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/DirichletBC.hpp>
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
22
#include <amdis/LocalAssemblerList.hpp>
23
#include <amdis/Marker.hpp>
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 grid
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
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
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
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
111 112
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath = {}, ColTreePath = {});
113

114
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
115 116
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath = {}, ColTreePath = {});
117
    /** @} */
118 119 120


    /// Adds an operator to \ref rhs.
121
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
122 123
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(Operator const& op, TreePath = {});
124

125
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
126 127
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath = {});
128
    /** @} */
129

130

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

137
  public:
138

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

144
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
145
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
146 147 148
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
149

150
    /// Writes output files.
151
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
152

Praetorius, Simon's avatar
Praetorius, Simon committed
153

154
  public: // get-methods
155

156 157 158
    /// Returns a pointer to system-matrix, \ref systemMatrix_
    std::shared_ptr<SystemMatrix> getSystemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
159

160
    /// Returns a pointer to the solution vector, \ref solution_
161
    std::shared_ptr<SystemVector> getSolutionVector() const
162
    {
163
      return solution_;
164 165 166 167 168 169 170
    }

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

174 175 176
    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
177
    {
178
      auto&& tp = makeTreePath(path);
179
      return makeDOFVectorView(*solution_, tp);
180
    }
181 182


183
    /// Return a point to the rhs system-vector, \ref rhs
184 185
    std::shared_ptr<SystemVector> getRhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
186 187


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

Praetorius, Simon's avatar
Praetorius, Simon committed
191
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
192
    {
193
      linearSolver_ = solver;
194 195
    }

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

Praetorius, Simon's avatar
Praetorius, Simon committed
199
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
200
    /// matrices and vectors, as well as the file-writer.
201
    void setGrid(Grid& grid)
202
    {
203 204
      grid_ = Dune::stackobject_to_shared_ptr(grid);

205
      createGlobalBasis();
206 207 208 209
      createMatricesAndVectors();
      createFileWriter();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
210 211
    /// Return the gridView of the leaf-level
    auto const& gridView() { return globalBasis_->gridView(); }
212

213
    /// Return the \ref feSpaces
214
    std::shared_ptr<GlobalBasis> const& globalBasis() { return globalBasis_; }
215 216


217
    /// Implementation of \ref ProblemStatBase::getName
218 219
    virtual std::string getName() const override
    {
220
      return name_;
221
    }
222

223
  protected: // initialization methods
224

225 226 227 228 229 230 231 232 233 234 235
    template <class T, class GV>
    using HasCreate = decltype(T::create(std::declval<GV>()));

    void createGlobalBasis()
    {
      createGlobalBasis(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
      initGlobalBasis(*globalBasis_);
    }

    void createGlobalBasis(std::true_type)
    {
236
      assert( bool(grid_) );
237
      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
238
      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
239 240 241 242 243 244 245
    }

    void createGlobalBasis(std::false_type)
    {
      error_exit("Can not create GlobalBasis from type. Pass a BasisCreator instead!");
    }

246
    void createGrid()
247
    {
248 249
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
250
      test_exit(!gridName_.empty(), "No grid 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 262
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
263 264 265 266
      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{});
267
    }
268

269 270
    void createMatricesAndVectors()
    {
271 272 273
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
274

275
      AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
276 277
      {
        std::string i = to_string(treePath);
278 279 280
        estimates_[i].resize(globalBasis_->gridView().indexSet().size(0));
        for (std::size_t j = 0; j < estimates_[i].size(); j++)
          estimates_[i][j] = 0.0; // TODO: Remove when estimate() is implemented
281 282 283
      });

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

286
    void createSolver()
287
    {
288
      std::string solverName = "cg";
289
      Parameters::get(name_ + "->solver->name", solverName);
290

291
      auto solverCreator
292
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
293

294
      linearSolver_ = solverCreator->create(name_ + "->solver");
295
    }
296

297
    void createMarker();
298

299
    void createFileWriter();
300

Praetorius, Simon's avatar
Praetorius, Simon committed
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
  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.
334
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
335 336


337
  private:
338
    /// Name of this problem.
339
    std::string name_;
340

341
    /// Grid of this problem.
342
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
343

344 345 346
    /// Number of grids
    int nGrids = 1;

347
    /// Name of the grid
348
    std::string gridName_ = "none";
349

350
    /// FE spaces of this problem.
351 352
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
353

354
    /// A FileWriter object
355
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
356

357
    /// Pointer to the adaptation markers
358
    std::list<std::shared_ptr<Marker<Traits> > > marker_;
359

360
    /// An object of the linearSolver Interface
361
    std::shared_ptr<LinearSolverType> linearSolver_;
362

363
    /// A block-matrix that is filled during assembling
364
    std::shared_ptr<SystemMatrix> systemMatrix_;
365

366
    /// A block-vector with the solution components
367
    std::shared_ptr<SystemVector> solution_;
368

369
    /// A vector with the local element error estimates
370 371
    /// for each node in the basis tree, indexed by [to_string(treePath)][element index]
    std::map<std::string, std::vector<double> > estimates_;
372

373
    /// A block-vector (load-vector) corresponding to the right.hand side
374
    /// of the equation, filled during assembling
375
    std::shared_ptr<SystemVector> rhs_;
376 377


378 379
  private: // some internal data-structures

380 381 382
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
383
  };
384

385 386 387 388 389 390 391 392 393 394 395 396 397 398 399
#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};
  }

400

401
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
402
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
403
#endif
404

405
} // end namespace AMDiS
406 407 408

#include "ProblemStat.inc.hpp"