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
#include <amdis/AdaptInfo.hpp>
#include <amdis/CreatorInterface.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
17
#include <amdis/CreatorMap.hpp>
18 19 20 21
#include <amdis/DirichletBC.hpp>
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
22
#include <amdis/LinearSolvers.hpp>
23
#include <amdis/LocalAssemblerList.hpp>
24
#include <amdis/Marker.hpp>
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/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 44 45 46
  // forward declaration
  template <class Traits>
  class ProblemInstat;

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

54 55
    friend class ProblemInstat<Traits>;

56
  public: // typedefs and static constants
57

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

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

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

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

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

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

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

103

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

114 115

    /// Adds an operator to \ref A.
116
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
117
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
118 119 120 121
    void addMatrixOperator(Operator const& op, RowTreePath row = {}, ColTreePath col = {})
    {
      systemMatrix_->addOperator(tag::element_operator<Element>{}, op, row, col);
    }
122

123
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
124
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
125 126 127 128
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath row = {}, ColTreePath col = {})
    {
      systemMatrix_->addOperator(tag::boundary_operator<typename GridView::Intersection>{b}, op, row, col);
    }
129
    /** @} */
130 131 132


    /// Adds an operator to \ref rhs.
133
    /** @{ */
Praetorius, Simon's avatar
Praetorius, Simon committed
134
    template <class Operator, class TreePath = RootTreePath>
135 136 137 138
    void addVectorOperator(Operator const& op, TreePath path = {})
    {
      rhs_->addOperator(tag::element_operator<Element>{}, op, path);
    }
139

140
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
141
    template <class Operator, class TreePath = RootTreePath>
142 143 144 145
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath path = {})
    {
      rhs_->addOperator(tag::boundary_operator<typename GridView::Intersection>{b}, op, path);
    }
146
    /** @} */
147

148

149
    /// Adds a Dirichlet boundary condition
150
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
151
    void addDirichletBC(Predicate const& predicate,
152
                        RowTreePath row, ColTreePath col,
Praetorius, Simon's avatar
Praetorius, Simon committed
153
                        Values const& values);
154

155
  public:
156

157
    /// Implementation of \ref ProblemStatBase::solve
158
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
159 160
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
161

162
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
163 164 165 166
    virtual void buildAfterAdapt(AdaptInfo& adaptInfo,
                                 Flag flag,
                                 bool asmMatrix = true,
                                 bool asmVector = true) override;
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 177 178 179 180 181 182 183 184
    /// 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_; }
185

186 187 188 189 190 191

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

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


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

208
    /// Set a new linear solver for the problem
Praetorius, Simon's avatar
Praetorius, Simon committed
209
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
210
    {
211
      linearSolver_ = solver;
212 213
    }

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

218
    /// Set the grid. Stores pointer and initializes feSpaces
219
    /// matrices and vectors, as well as the file-writer.
220
    void setGrid(std::shared_ptr<Grid> const& grid)
221
    {
222
      grid_ = grid;
223

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

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

232 233 234 235 236
    /// Return the \ref globalBasis_
    GlobalBasis const& globalBasis() const { return *globalBasis_; }

    /// Return the \ref globalBasis_
    GlobalBasis& globalBasis() { 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 248 249 250 251 252 253 254 255 256 257
    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)
    {
258
      assert( bool(grid_) );
259
      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
260
      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
261 262 263 264
    }

    void createGlobalBasis(std::false_type)
    {
265
      error_exit("Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
266 267
    }

268
    void createGrid()
269
    {
270 271
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
272
      test_exit(!gridName_.empty(), "No grid name specified for '{}->mesh'!", name_);
273

274
      grid_ = MeshCreator<Grid>::create(gridName_);
275

276
      msg("Create grid:");
277 278 279
      msg("#elements = {}"   , grid_->size(0));
      msg("#faces/edges = {}", grid_->size(1));
      msg("#vertices = {}"   , grid_->size(dim));
280
      msg("");
281
    }
282

283 284
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
285
      constraints_.init(*globalBasis_, *globalBasis_);
286
    }
287

288 289
    void createMatricesAndVectors()
    {
290 291 292
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_);
      solution_ = std::make_shared<SystemVector>(*globalBasis_);
      rhs_ = std::make_shared<SystemVector>(*globalBasis_);
293

294 295
      auto localView = globalBasis_->localView();
      AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
296 297
      {
        std::string i = to_string(treePath);
298 299 300
        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
301
      });
302
    }
303

304
    void createSolver()
305
    {
306
      std::string solverName = "cg";
307
      Parameters::get(name_ + "->solver->name", solverName);
308

309
      auto solverCreator
310
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
311

312
      linearSolver_ = solverCreator->create(name_ + "->solver");
313
    }
314

315
    void createMarker();
316

317
    void createFileWriter();
318

Praetorius, Simon's avatar
Praetorius, Simon committed
319

320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339
  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.
340
    virtual Flag adaptGrid(AdaptInfo& adaptInfo) override;
341 342

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


346
  private:
347
    /// Name of this problem.
348
    std::string name_;
349

350
    /// Grid of this problem.
351
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
352

353 354 355
    /// Number of grids
    int nGrids = 1;

356
    /// Name of the grid
357
    std::string gridName_ = "none";
358

359
    /// FE spaces of this problem.
360
    std::shared_ptr<GlobalBasis> globalBasis_;
361

362
    /// A FileWriter object
363
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
364

365
    /// Pointer to the adaptation markers
366
    std::list<std::shared_ptr<Marker<Traits> > > marker_;
367

368
    /// An object of the linearSolver Interface
369
    std::shared_ptr<LinearSolverType> linearSolver_;
370

371
    /// A block-matrix that is filled during assembling
372
    std::shared_ptr<SystemMatrix> systemMatrix_;
373

374
    /// A block-vector with the solution components
375
    std::shared_ptr<SystemVector> solution_;
376

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

381
    /// A block-vector (load-vector) corresponding to the right.hand side
382
    /// of the equation, filled during assembling
383
    std::shared_ptr<SystemVector> rhs_;
384 385


386 387
  private: // some internal data-structures

388
    Constraints<GlobalBasis, GlobalBasis> constraints_;
389
  };
390

391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
#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};
  }

406

407
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
Praetorius, Simon's avatar
Praetorius, Simon committed
408 409 410
  // extern template class ProblemStat<YaspGridBasis<2,1>>;
  // extern template class ProblemStat<YaspGridBasis<2,2>>;
  // extern template class ProblemStat<YaspGridBasis<2,1,2>>;
411
#endif
412

413
} // end namespace AMDiS
414 415 416

#include "ProblemStat.inc.hpp"