ProblemStat.hpp 12.7 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/Marker.hpp>
23 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/TypeDefs.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
  template <class Traits>
43
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
44
      : public ProblemStatBase
45
      , public StandardProblemIteration
46
  {
47
    using Self = ProblemStat;
48

49
  public: // typedefs and static constants
50

51
    using GlobalBasis = typename Traits::GlobalBasis;
52 53 54
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
    using Element  = typename GridView::template Codim<0>::Entity;
55
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
56

57
    /// Dimension of the grid
58
    static constexpr int dim = Grid::dimension;
59

60
    /// Dimension of the world
61
    static constexpr int dow = Grid::dimensionworld;
62

63
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
64
    using SystemVector = DOFVector<Traits, double>;
65

66
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
67

68
  public:
69 70
    /**
     * \brief Constructor. Takes the name of the problem that is used to
71
     * access values corresponding to this problem in the parameter file.
72
     **/
73
    explicit ProblemStat(std::string name)
74
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
75
      , name_(std::move(name))
76
    {}
77

78 79
    /// 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
80
    ProblemStat(std::string name, Grid& grid)
81
      : ProblemStat(std::move(name))
82
    {
83 84 85 86 87 88 89 90 91 92 93
      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_);
94 95
    }

96

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

107 108

    /// Adds an operator to \ref A.
109
    /** @{ */
110 111
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath, ColTreePath,
112
                           double* factor = nullptr, double* estFactor = nullptr);
113

114
    // operator evaluated on the boundary
115 116
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath, ColTreePath,
117
                           double* factor = nullptr, double* estFactor = nullptr);
118
    /** @} */
119 120 121


    /// Adds an operator to \ref rhs.
122
    /** @{ */
123 124
    template <class Operator, class TreePath>
    void addVectorOperator(Operator const& op, TreePath,
125
                           double* factor = nullptr, double* estFactor = nullptr);
126

127
    // operator evaluated on the boundary
128 129
    template <class Operator, class TreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath,
130
                           double* factor = nullptr, double* estFactor = nullptr);
131
    /** @} */
132

133

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

140
  public:
141

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

147
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
148
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
149 150 151
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
152

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

Praetorius, Simon's avatar
Praetorius, Simon committed
156

157
  public: // get-methods
158

159 160 161
    /// Returns a pointer to system-matrix, \ref systemMatrix_
    std::shared_ptr<SystemMatrix> getSystemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
162

163
    /// Returns a pointer to the solution vector, \ref solution_
164
    std::shared_ptr<SystemVector> getSolutionVector() const
165
    {
166
      return solution_;
167 168 169 170 171 172 173
    }

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

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


186
    /// Return a point to the rhs system-vector, \ref rhs
187 188
    std::shared_ptr<SystemVector> getRhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
189 190


191
    /// Return a pointer to the linear solver, \ref linearSolver
192
    std::shared_ptr<LinearSolverType> getSolver() { return linearSolver_; }
193

Praetorius, Simon's avatar
Praetorius, Simon committed
194
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
195
    {
196
      linearSolver_ = solver;
197 198
    }

199
    /// Return a pointer to the grid, \ref grid
200
    std::shared_ptr<Grid> getGrid() { return grid_; }
201

202
    /// Set the grid. Stores pointer to passed reference and initializes feSpaces
203
    /// matrices and vectors, as well as the file-writer.
204
    void setGrid(Grid& grid)
205
    {
206 207
      grid_ = Dune::stackobject_to_shared_ptr(grid);

208
      createGlobalBasis();
209 210 211 212
      createMatricesAndVectors();
      createFileWriter();
    }

213

214
    /// Return the \ref feSpaces
215
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
216 217


218
    /// Implementation of \ref ProblemStatBase::getName
219 220
    virtual std::string getName() const override
    {
221
      return name_;
222
    }
223

224
  protected: // initialization methods
225

226
    void createGrid()
227
    {
228 229
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
230
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
231

232
      grid_ = MeshCreator<Grid>::create(gridName_);
233

234
      msg("Create grid:");
235 236 237
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
238
      msg("");
239
    }
240

241
    void createGlobalBasis()
242
    {
243 244
      globalBasis_ = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid_->leafGridView()));
      initGlobalBasis(*globalBasis_);
245 246 247 248
    }

    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
249 250 251 252
      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{});
253
    }
254

255 256
    void createMatricesAndVectors()
    {
257 258 259
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
260

261
      AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
262 263
      {
        std::string i = to_string(treePath);
264 265 266
        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
267 268 269
      });

      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
270
    }
271

272
    void createSolver()
273
    {
274
      std::string solverName = "cg";
275
      Parameters::get(name_ + "->solver->name", solverName);
276

277
      auto solverCreator
278
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
279

280
      linearSolver_ = solverCreator->create(name_ + "->solver");
281
    }
282

283
    void createMarker();
284

285
    void createFileWriter();
286

Praetorius, Simon's avatar
Praetorius, Simon committed
287

288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319
  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::buildBeforeRefine.
    virtual void buildBeforeRefine(AdaptInfo&, Flag) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::buildBeforeCoarsen.
    virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::estimate.
    virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::refineMesh.
    virtual Flag refineMesh(AdaptInfo& adaptInfo) override { return 0; }

    /// Implementation of \ref ProblemStatBase::coarsenMesh.
    virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return 0; }

    /// Implementation of \ref ProblemStatBase::markElements.
320
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
321 322


323
  private:
324
    /// Name of this problem.
325
    std::string name_;
326

327
    /// Grid of this problem.
328
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
329

330 331 332
    /// Number of grids
    int nGrids = 1;

333
    /// Name of the grid
334
    std::string gridName_ = "none";
335

336
    /// FE spaces of this problem.
337 338
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
339

340
    /// A FileWriter object
341
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
342

343
    /// Pointer to the adaptation markers
344
    std::list<std::shared_ptr<Marker<Traits> > > marker_;
345

346
    /// An object of the linearSolver Interface
347
    std::shared_ptr<LinearSolverType> linearSolver_;
348

349
    /// A block-matrix that is filled during assembling
350
    std::shared_ptr<SystemMatrix> systemMatrix_;
351

352
    /// A block-vector with the solution components
353
    std::shared_ptr<SystemVector> solution_;
354

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

359
    /// A block-vector (load-vector) corresponding to the right.hand side
360
    /// of the equation, filled during assembling
361
    std::shared_ptr<SystemVector> rhs_;
362 363


364 365
  private: // some internal data-structures

366 367 368
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
369
  };
370

371 372 373 374 375 376 377 378 379 380 381 382 383 384 385
#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};
  }

386

387
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
388
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
389
#endif
390

391
} // end namespace AMDiS
392 393 394

#include "ProblemStat.inc.hpp"