ProblemStat.hpp 13.5 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
#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>
35
#include <amdis/gridfunctions/DiscreteFunction.hpp>
36 37 38 39
#include <amdis/gridfunctions/DOFVectorView.hpp>

#include <amdis/io/FileWriterInterface.hpp>

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

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

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

56 57
    friend class ProblemInstat<Traits>;

58
  public: // typedefs and static constants
59

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

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

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

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

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

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

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

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
98 99
    ProblemStat(std::string const& name, Grid& grid, GlobalBasis& globalBasis)
      : ProblemStat(name, grid)
100 101 102
    {
      globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
      initGlobalBasis(*globalBasis_);
103
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
104

105

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

116 117

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

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


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

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

150

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

157
  public:
158

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
173

174
  public: // get-methods
175

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

188 189 190 191 192 193

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

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


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

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
231 232 233 234 235 236 237 238 239 240 241 242
    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
243
    /// Return the gridView of the leaf-level
244
    GridView const& gridView() { return globalBasis_->gridView(); }
245

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

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


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

259
  protected: // initialization methods
260

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

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

282
    void createGrid()
283
    {
284 285
      Parameters::get(name_ + "->mesh", gridName_);
      grid_ = MeshCreator<Grid>::create(gridName_);
286

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

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

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

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

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

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

323
      linearSolver_ = solverCreator->create(name_ + "->solver");
324
    }
325

326
    void createMarker();
327

328
    void createFileWriter();
329

Praetorius, Simon's avatar
Praetorius, Simon committed
330

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

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


357
  private:
358
    /// Name of this problem.
359
    std::string name_;
360

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

364 365 366
    /// Number of grids
    int nGrids = 1;

367
    /// Name of the grid
368
    std::string gridName_ = "mesh";
369

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

373
    /// A FileWriter object
374
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
375

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

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

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

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

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

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

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


400 401
  private: // some internal data-structures

402
    Constraints<GlobalBasis, GlobalBasis> constraints_;
403
  };
404

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

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

420

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

427
} // end namespace AMDiS
428 429 430

#include "ProblemStat.inc.hpp"