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

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

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

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

53 54
    friend class ProblemInstat<Traits>;

55
  public: // typedefs and static constants
56

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

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

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

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

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

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

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

102

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

113 114

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

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


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

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

135

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

142
  public:
143

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
158

159
  public: // get-methods
160

161 162 163 164 165 166 167 168 169 170 171
    /// 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_; }
172

173 174 175 176 177 178

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

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


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

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

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

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

211
      createGlobalBasis();
212 213 214 215
      createMatricesAndVectors();
      createFileWriter();
    }

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

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

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


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

232
  protected: // initialization methods
233

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

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

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

261
      grid_ = MeshCreator<Grid>::create(gridName_);
262

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

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

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

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

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

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

301
      linearSolver_ = solverCreator->create(name_ + "->solver");
302
    }
303

304
    void createMarker();
305

306
    void createFileWriter();
307

Praetorius, Simon's avatar
Praetorius, Simon committed
308

309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328
  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.
329
    virtual Flag adaptGrid(AdaptInfo& adaptInfo) override;
330 331

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


335
  private:
336
    /// Name of this problem.
337
    std::string name_;
338

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

342 343 344
    /// Number of grids
    int nGrids = 1;

345
    /// Name of the grid
346
    std::string gridName_ = "none";
347

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

352
    /// A FileWriter object
353
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
354

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

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

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

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

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

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


376 377
  private: // some internal data-structures

378 379 380
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
381
  };
382

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

398

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

403
} // end namespace AMDiS
404 405 406

#include "ProblemStat.inc.hpp"