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
#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/CreatorInterface.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
18
#include <amdis/CreatorMap.hpp>
19 20 21 22
#include <amdis/DirichletBC.hpp>
#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 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/Utility.hpp>

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

#include <amdis/io/FileWriterInterface.hpp>

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

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

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

55 56
    friend class ProblemInstat<Traits>;

57
  public: // typedefs and static constants
58

59
    using GlobalBasis = typename Traits::GlobalBasis;
60 61 62
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
    using Element  = typename GridView::template Codim<0>::Entity;
63
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
64

65
    /// Dimension of the grid
66
    static constexpr int dim = Grid::dimension;
67

68
    /// Dimension of the world
69
    static constexpr int dow = Grid::dimensionworld;
70

71
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
72
    using SystemVector = DOFVector<GlobalBasis, double>;
73

74
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
75

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

86 87
    /// 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
88
    ProblemStat(std::string name, Grid& grid)
89
      : ProblemStat(std::move(name))
90
    {
91 92 93 94 95 96 97 98 99 100 101
      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_);
102
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
103

104

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

115 116

    /// Adds an operator to \ref A.
117
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
118 119
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath = {}, ColTreePath = {});
120

121
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
122 123
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath = {}, ColTreePath = {});
124
    /** @} */
125 126 127


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

132
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
133 134
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath = {});
135
    /** @} */
136

137

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

144
  public:
145

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

151
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
152 153 154 155
    virtual void buildAfterAdapt(AdaptInfo& adaptInfo,
                                 Flag flag,
                                 bool asmMatrix = true,
                                 bool asmVector = true) override;
156

157
    /// Writes output files.
158
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
159

Praetorius, Simon's avatar
Praetorius, Simon committed
160

161
  public: // get-methods
162

163 164 165 166 167 168 169 170 171 172 173
    /// Returns a reference to system-matrix, \ref systemMatrix_
    SystemMatrix&       getSystemMatrix()       { return *systemMatrix_; }
    SystemMatrix const& getSystemMatrix() const { return *systemMatrix_; }

    /// Returns a reference to the solution vector, \ref solution_
    SystemVector&       getSolutionVector()       { return *solution_; }
    SystemVector const& getSolutionVector() const { return *solution_; }

    /// Return a reference to the rhs system-vector, \ref rhs
    SystemVector&       getRhsVector()       { return *rhs_; }
    SystemVector const& getRhsVector() const { return *rhs_; }
174

175 176 177 178 179 180

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

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


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

197
    /// Set a new linear solver for the problem
Praetorius, Simon's avatar
Praetorius, Simon committed
198
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
199
    {
200
      linearSolver_ = solver;
201 202
    }

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

207
    /// Set the grid. Stores pointer and initializes feSpaces
208
    /// matrices and vectors, as well as the file-writer.
209
    void setGrid(std::shared_ptr<Grid> const& grid)
210
    {
211
      grid_ = grid;
212

213
      createGlobalBasis();
214 215 216 217
      createMatricesAndVectors();
      createFileWriter();
    }

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

221 222 223 224 225
    /// Return the \ref globalBasis_
    GlobalBasis const& globalBasis() const { return *globalBasis_; }

    /// Return the \ref globalBasis_
    GlobalBasis& globalBasis() { return *globalBasis_; }
226 227


228
    /// Implementation of \ref ProblemStatBase::getName
229 230
    virtual std::string getName() const override
    {
231
      return name_;
232
    }
233

234
  protected: // initialization methods
235

236 237 238 239 240 241 242 243 244 245 246
    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)
    {
247
      assert( bool(grid_) );
248
      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
249
      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
250 251 252 253
    }

    void createGlobalBasis(std::false_type)
    {
254
      error_exit("Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
255 256
    }

257
    void createGrid()
258
    {
259 260
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
261
      test_exit(!gridName_.empty(), "No grid name specified for '{}->mesh'!", name_);
262

263
      grid_ = MeshCreator<Grid>::create(gridName_);
264

265
      msg("Create grid:");
266 267 268
      msg("#elements = {}"   , grid_->size(0));
      msg("#faces/edges = {}", grid_->size(1));
      msg("#vertices = {}"   , grid_->size(dim));
269
      msg("");
270
    }
271

272 273
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
274 275 276 277
      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{});
278
    }
279

280 281
    void createMatricesAndVectors()
    {
282 283 284
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_);
      solution_ = std::make_shared<SystemVector>(*globalBasis_);
      rhs_ = std::make_shared<SystemVector>(*globalBasis_);
285

286
      AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
287 288
      {
        std::string i = to_string(treePath);
289 290 291
        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
292
      });
293
    }
294

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

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

303
      linearSolver_ = solverCreator->create(name_ + "->solver");
304
    }
305

306
    void createMarker();
307

308
    void createFileWriter();
309

Praetorius, Simon's avatar
Praetorius, Simon committed
310

311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330
  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::estimate.
    virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }

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

    /// 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
Praetorius, Simon's avatar
Praetorius, Simon committed
402 403 404
  // extern template class ProblemStat<YaspGridBasis<2,1>>;
  // extern template class ProblemStat<YaspGridBasis<2,2>>;
  // extern template class ProblemStat<YaspGridBasis<2,1,2>>;
405
#endif
406

407
} // end namespace AMDiS
408 409 410

#include "ProblemStat.inc.hpp"