ProblemStat.hpp 13.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
#include <amdis/AdaptInfo.hpp>
#include <amdis/CreatorInterface.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
17
#include <amdis/CreatorMap.hpp>
18
#include <amdis/DirichletBC.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
19
//#include <amdis/Estimator.hpp>
20 21 22
#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 const& name)
82
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
83
      , name_(name)
84
    {}
85

86 87
    /// Constructor taking additionally a reference to a grid that is used
    /// instead of the default created grid, \ref ProblemStat
88 89
    ProblemStat(std::string const& name, Grid& grid)
      : ProblemStat(name)
90
    {
91
      grid_ = Dune::stackobject_to_shared_ptr(grid);
92
      Parameters::get(name + "->mesh", gridName_);
93 94 95 96
    }

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
97 98
    ProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
      : ProblemStat(name, grid)
99 100 101
    {
      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
    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
119 120 121 122
    void addMatrixOperator(Operator const& op, RowTreePath row = {}, ColTreePath col = {})
    {
      systemMatrix_->addOperator(tag::element_operator<Element>{}, op, row, col);
    }
123

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


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

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

149

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

156
  public:
157

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

163
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
164 165 166 167
    virtual void buildAfterAdapt(AdaptInfo& adaptInfo,
                                 Flag flag,
                                 bool asmMatrix = true,
                                 bool asmVector = true) override;
168

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

Praetorius, Simon's avatar
Praetorius, Simon committed
172

173
  public: // get-methods
174

175 176 177 178 179 180 181 182 183 184 185
    /// 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_; }
186

187 188 189 190 191 192

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

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


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

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

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

219
    /// Set the grid. Stores pointer and initializes feSpaces
Praetorius, Simon's avatar
Praetorius, Simon committed
220
    /// matrices and vectors, as well as markers and file-writers.
221
    void setGrid(std::shared_ptr<Grid> const& grid)
222
    {
223
      grid_ = grid;
224
      createGlobalBasis();
225
      createMatricesAndVectors();
Praetorius, Simon's avatar
Praetorius, Simon committed
226
      createMarker();
227 228 229
      createFileWriter();
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
230 231 232 233 234 235 236 237 238 239 240 241
    void addMarker(std::shared_ptr<Marker<Grid>> const& marker)
    {
      marker_.push_back(marker);
      if (marker_.size() > 1)
        marker_.back()->setMaximumMarking(true);
    }

    void addMarker(Marker<Grid>& marker)
    {
      addMarker(Dune::stackobject_to_shared_ptr(marker));
    }

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

245 246 247 248 249
    /// Return the \ref globalBasis_
    GlobalBasis const& globalBasis() const { return *globalBasis_; }

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


252
    /// Implementation of \ref ProblemStatBase::getName
253
    virtual std::string const& getName() const override
254
    {
255
      return name_;
256
    }
257

258
  protected: // initialization methods
259

260 261 262 263 264 265 266 267 268 269 270
    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)
    {
271
      assert( bool(grid_) );
272
      static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
273
      globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
274 275 276 277
    }

    void createGlobalBasis(std::false_type)
    {
278
      error_exit("Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
279 280
    }

281
    void createGrid()
282
    {
283 284
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
285
      test_exit(!gridName_.empty(), "No grid name specified for '{}->mesh'!", name_);
286

287
      grid_ = MeshCreator<Grid>::create(gridName_);
288

289
      msg("Create grid:");
290 291 292
      msg("#elements = {}"   , grid_->size(0));
      msg("#faces/edges = {}", grid_->size(1));
      msg("#vertices = {}"   , grid_->size(dim));
293
      msg("");
294
    }
295

296 297
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
298
      constraints_.init(*globalBasis_, *globalBasis_);
299
    }
300

301 302
    void createMatricesAndVectors()
    {
303 304 305
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_);
      solution_ = std::make_shared<SystemVector>(*globalBasis_);
      rhs_ = std::make_shared<SystemVector>(*globalBasis_);
306

307 308
      auto localView = globalBasis_->localView();
      AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
309 310
      {
        std::string i = to_string(treePath);
311 312 313
        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
314
      });
315
    }
316

317
    void createSolver()
318
    {
319
      std::string solverName = "cg";
320
      Parameters::get(name_ + "->solver->name", solverName);
321

322
      auto solverCreator
323
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
324

325
      linearSolver_ = solverCreator->create(name_ + "->solver");
326
    }
327

328
    void createMarker();
329

330
    void createFileWriter();
331

Praetorius, Simon's avatar
Praetorius, Simon committed
332

333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352
  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.
353
    virtual Flag adaptGrid(AdaptInfo& adaptInfo) override;
354 355

    /// Implementation of \ref ProblemStatBase::markElements.
356
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
357 358


359
  private:
360
    /// Name of this problem.
361
    std::string name_;
362

363
    /// Grid of this problem.
364
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
365

366 367 368
    /// Number of grids
    int nGrids = 1;

369
    /// Name of the grid
370
    std::string gridName_ = "none";
371

372
    /// FE spaces of this problem.
373
    std::shared_ptr<GlobalBasis> globalBasis_;
374

375
    /// A FileWriter object
376
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
377

378
    /// Pointer to the adaptation markers
Praetorius, Simon's avatar
Praetorius, Simon committed
379 380 381 382
    std::list<std::shared_ptr<Marker<Grid> > > marker_;

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

384
    /// An object of the linearSolver Interface
385
    std::shared_ptr<LinearSolverType> linearSolver_;
386

387
    /// A block-matrix that is filled during assembling
388
    std::shared_ptr<SystemMatrix> systemMatrix_;
389

390
    /// A block-vector with the solution components
391
    std::shared_ptr<SystemVector> solution_;
392

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

397
    /// A block-vector (load-vector) corresponding to the right.hand side
398
    /// of the equation, filled during assembling
399
    std::shared_ptr<SystemVector> rhs_;
400 401


402 403
  private: // some internal data-structures

404
    Constraints<GlobalBasis, GlobalBasis> constraints_;
405
  };
406

407 408 409
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
  // Deduction rule
  template <class Grid, class GlobalBasis>
410
  ProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
411 412 413 414 415 416
    -> ProblemStat<DefaultProblemTraits<GlobalBasis>>;
#endif

  // Generator for ProblemStat with given Grid and GlobalBasis
  template <class Grid, class GlobalBasis>
  ProblemStat<DefaultProblemTraits<GlobalBasis>>
417
  makeProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
418
  {
419
    return {name, grid, globalBasis};
420 421
  }

422

423
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
Praetorius, Simon's avatar
Praetorius, Simon committed
424 425 426
  // extern template class ProblemStat<YaspGridBasis<2,1>>;
  // extern template class ProblemStat<YaspGridBasis<2,2>>;
  // extern template class ProblemStat<YaspGridBasis<2,1,2>>;
427
#endif
428

429
} // end namespace AMDiS
430 431 432

#include "ProblemStat.inc.hpp"