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

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

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

50
  public: // typedefs and static constants
51

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

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

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

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

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

69
  public:
Praetorius, Simon's avatar
Praetorius, Simon committed
70
#if 0
71 72
    /**
     * \brief Constructor. Takes the name of the problem that is used to
73
     * access values corresponding to this problem in the parameter file.
74
     **/
75
    explicit ProblemStat(std::string name)
76
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
Praetorius, Simon's avatar
Praetorius, Simon committed
77
      , name(std::move(name))
78
    {}
79

80 81
    /// 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
82
    ProblemStat(std::string name, Grid& grid)
83
      : ProblemStat(std::move(name))
84
    {
85 86 87 88 89 90 91 92 93 94 95
      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_);
96
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
97 98 99 100 101 102 103 104 105 106 107 108 109 110
#endif

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
    ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      , name(std::move(name))
    {
      gridName = "";
      Parameters::get(name + "->mesh", gridName);

      globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
      initGlobalBasis(*globalBasis_);
    }
111

112

113
    /**
114
     * \brief Initialisation of the problem.
115
     *
116 117 118
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
119
    void initialize(Flag initFlag,
120
                    Self* adoptProblem = nullptr,
Praetorius, Simon's avatar
Praetorius, Simon committed
121
                    Flag adoptFlag = INIT_NOTHING);
122

123 124

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

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


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

140
    // operator evaluated on the boundary
Praetorius, Simon's avatar
Praetorius, Simon committed
141 142
    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath = {});
143
    /** @} */
144

145

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

152
  public:
153

154
    /// Implementation of \ref ProblemStatBase::solve
155
    virtual void solve(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
156 157
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;
158

159
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
160
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
Praetorius, Simon's avatar
Praetorius, Simon committed
161 162 163
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;
164

165
    /// Writes output files.
166
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);
167

Praetorius, Simon's avatar
Praetorius, Simon committed
168

169
  public: // get-methods
170

171 172 173
    /// Returns a pointer to system-matrix, \ref systemMatrix_
    std::shared_ptr<SystemMatrix> getSystemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
174

175
    /// Returns a pointer to the solution vector, \ref solution_
176
    std::shared_ptr<SystemVector> getSolutionVector() const
177
    {
178
      return solution_;
179 180 181 182 183 184 185
    }

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

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


198
    /// Return a point to the rhs system-vector, \ref rhs
199 200
    std::shared_ptr<SystemVector> getRhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
201 202


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

Praetorius, Simon's avatar
Praetorius, Simon committed
206
    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
207
    {
208
      linearSolver_ = solver;
209 210
    }

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

Praetorius, Simon's avatar
Praetorius, Simon committed
214 215
#if 0
    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
216
    /// matrices and vectors, as well as the file-writer.
217
    void setGrid(Grid& grid)
218
    {
219 220
      grid_ = Dune::stackobject_to_shared_ptr(grid);

221
      createGlobalBasis();
222 223 224
      createMatricesAndVectors();
      createFileWriter();
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
225
#endif
226

Praetorius, Simon's avatar
Praetorius, Simon committed
227 228
    /// Return the gridView of the leaf-level
    auto const& gridView() { return globalBasis_->gridView(); }
229

230
    /// Return the \ref feSpaces
231
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
232 233


234
    /// Implementation of \ref ProblemStatBase::getName
235 236
    virtual std::string getName() const override
    {
237
      return name_;
238
    }
239

240
  protected: // initialization methods
241

242
    void createGrid()
243
    {
244 245
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
246
      test_exit(!gridName_.empty(), "No grid name specified for '", name_, "->mesh'!");
247

248
      grid_ = MeshCreator<Grid>::create(gridName_);
249

250
      msg("Create grid:");
251 252 253
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
254
      msg("");
255
    }
256

257 258
    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
259
      localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
Praetorius, Simon's avatar
Praetorius, Simon committed
260 261 262
      matrixOperators.init(localView_->tree(), tag::store{});
      rhsOperators.init(localView_->tree(), tag::store{});
      constraints.init(localView_->tree(), tag::store{});
263
    }
264

265 266
    void createMatricesAndVectors()
    {
267 268 269
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");
270

271
      AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
272 273
      {
        std::string i = to_string(treePath);
274 275 276
        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
277 278 279
      });

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

282
    void createSolver()
283
    {
284
      std::string solverName = "cg";
285
      Parameters::get(name_ + "->solver->name", solverName);
286

287
      auto solverCreator
288
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));
289

290
      linearSolver_ = solverCreator->create(name_ + "->solver");
291
    }
292

293
    void createMarker();
294

295
    void createFileWriter();
296

Praetorius, Simon's avatar
Praetorius, Simon committed
297

298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329
  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.
330
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
331 332


333
  private:
334
    /// Name of this problem.
335
    std::string name_;
336

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

340 341 342
    /// Number of grids
    int nGrids = 1;

343
    /// Name of the grid
344
    std::string gridName_ = "none";
345

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

350
    /// A FileWriter object
351
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;
352

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

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

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

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

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

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


374 375
  private: // some internal data-structures

376 377 378
    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
379
  };
380

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

396

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

401
} // end namespace AMDiS
402 403 404

#include "ProblemStat.inc.hpp"