ProblemStat.hpp 12.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 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
#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>
#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>

#include <amdis/utility/TreePath.hpp>
37

38
namespace AMDiS
39
{
40
  template <class Traits>
41
  class ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
42
      : public ProblemStatBase
43
      , public StandardProblemIteration
44
  {
45
    using Self = ProblemStat;
46

47
  public: // typedefs and static constants
48

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

    /// Dimension of the mesh
56
    static constexpr int dim = Grid::dimension;
57

58
    /// Dimension of the world
59
    static constexpr int dow = Grid::dimensionworld;
60

61
    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
62
    using SystemVector = DOFVector<Traits, double>;
63

64
    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
65

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

76 77
    /// Constructor taking additionally a reference to a mesh that is used
    /// instead of the default created mesh, \ref ProblemStat
Praetorius, Simon's avatar
Praetorius, Simon committed
78
    ProblemStat(std::string name, Grid& grid)
79
      : ProblemStat(std::move(name))
80
    {
81 82 83 84 85 86 87 88 89 90 91
      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_);
92 93
    }

94

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

105 106

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

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


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

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

131

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

138
  public:
139

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


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

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 203
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// 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 213
      createMatricesAndVectors();
      createFileWriter();
    }


214
    /// Return the gridView of the leaf-level
215 216 217 218
    auto leafGridView() { return grid_->leafGridView(); }

    /// Return the gridView of levle `level`
    auto levelGridView(int level) { return grid_->levelGridView(level); }
219

220
    /// Return the \ref feSpaces
221
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
222 223


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

230
  protected: // initialization methods
231

232
    void createGrid()
233
    {
234 235 236
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
      test_exit(!gridName_.empty(), "No mesh name specified for '", name_, "->mesh'!");
237

238
      grid_ = MeshCreator<Grid>::create(gridName_);
239

240
      msg("Create grid:");
241 242 243
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
244
      msg("");
245
    }
246

247
    void createGlobalBasis()
248
    {
249 250
      globalBasis_ = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid_->leafGridView()));
      initGlobalBasis(*globalBasis_);
251 252 253 254
    }

    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
255 256 257 258
      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{});
259
    }
260

261 262
    void createMatricesAndVectors()
    {
263 264 265
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
266
    }
267

268
    void createSolver()
269
    {
270
      std::string solverName = "cg";
271
      Parameters::get(name_ + "->solver->name", solverName);
272

273
      auto solverCreator
274
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
275

276
      linearSolver_ = solverCreator->create(name_ + "->solver");
277
    }
278

279
    void createFileWriter();
280

Praetorius, Simon's avatar
Praetorius, Simon committed
281

282 283 284 285 286 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
  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.
    virtual Flag markElements(AdaptInfo& adaptInfo) override { return 0; }


317
  private:
318
    /// Name of this problem.
319
    std::string name_;
320

321
    /// Grid of this problem.
322
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems
323

324
    /// Name of the mesh
325
    std::string gridName_ = "none";
326

327
    /// FE spaces of this problem.
328 329
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;
330

331
    /// A FileWriter object
332
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
333

334
    /// An object of the linearSolver Interface
335
    std::shared_ptr<LinearSolverType> linearSolver_;
336

337
    /// A block-matrix that is filled during assembling
338
    std::shared_ptr<SystemMatrix> systemMatrix_;
339

340
    /// A block-vector with the solution components
341
    std::shared_ptr<SystemVector> solution_;
342 343

    /// A block-vector (load-vector) corresponding to the right.hand side
344
    /// of the equation, filled during assembling
345
    std::shared_ptr<SystemVector> rhs_;
346 347


348 349
  private: // some internal data-structures

350 351 352
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
353
  };
354

355 356 357 358 359 360 361 362 363 364 365 366 367 368 369
#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};
  }

370

371
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
372
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
373
#endif
374

375
} // end namespace AMDiS
376 377 378

#include "ProblemStat.inc.hpp"