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

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

97

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

108 109

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

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


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

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

134

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

141
  public:
142

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
157

158
  public: // get-methods
159

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

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

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

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


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


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

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

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

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

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

214

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


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

225
  protected: // initialization methods
226

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

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

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

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

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

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

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

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

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

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

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

284
    void createMarker();
285

286
    void createFileWriter();
287

Praetorius, Simon's avatar
Praetorius, Simon committed
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 320
  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.
321
    virtual Flag markElements(AdaptInfo& adaptInfo) override;
322 323


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

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

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

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

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

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

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

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

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

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

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

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


365 366
  private: // some internal data-structures

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

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

387

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

392
} // end namespace AMDiS
393 394 395

#include "ProblemStat.inc.hpp"